!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

!=======================================================================
! FVCOM Sediment Module 
!
! Copyright:    2005(c)
!
! THIS IS A DEMONSTRATION RELEASE. THE AUTHOR(S) MAKE NO REPRESENTATION
! ABOUT THE SUITABILITY OF THIS SOFTWARE FOR ANY OTHER PURPOSE. IT IS
! PROVIDED "AS IS" WITHOUT EXPRESSED OR IMPLIED WARRANTY.
!
! THIS ORIGINAL HEADER MUST BE MAINTAINED IN ALL DISTRIBUTED
! VERSIONS.
!
! Contact:      G. Cowles 
!               School for Marine Science and Technology, Umass-Dartmouth
!
! Based on the Community Sediment Transport Model (CSTM) as implemented
!     in ROMS by J. Warner (USGS)
!
! Comments:     Sediment Dynamics Module 
!
! Current FVCOM dependency
!
!   init_sed.F:  - user defined sediment model initial conditions
!   mod_ncdio.F: - netcdf output includes concentration/bottom/bed fields
!   fvcom.F:     - main calls sediment setup 
!   internal.F:  - calls sediment advance
!
! History
!   Feb 7, 2008: added initialization of bottom(:,:) to 0 (w/ T. Hamada)
!              : fixed loop bounds in hot start and archive for conc (w/ T. Hamada)
!              : added comments describing theoretical bases of dynamics
!   Feb 14,2008: added non-constant settling velocity for cohesive sediments (w/ T. Hamada) 
!              : updated vertical flux routine to handle non-constant vertical velocity (w/ T. Hamada)
!              : added a user-defined routine to calculate settling velocity based on concentration (w/ T. Hamada)
!              : added a user-defined routine to calculate erosion for a general case (w/ T. Hamada)
!
!  PLEASE NOTE!!!!!!!!!!! 
!  Do NOT USE INTEL FORTRAN COMPILER VERSION 11.0 IT HAS KNOWN BUGS WHEN DEALING WITH TYPES WITH ALLOCATABLE
!    COMPONENTS.  YOU WILL SEE WEIRD BEHAVIOR.  VERSION 11.1 IS OK.
!   
!
!  Later
!   1.) Modify vertical flux routines to work with general vertical coordinate
!   2.) Add divergence term for bedload transport calc 
!   3.) Add ripple roughness calculation
!   4.) Add morphological change (bathymetry + vertical velocity condition) 
!   5.) Eliminate excess divisions and recalcs
!
!=======================================================================
Module Mod_Sed  
#if defined (SEDIMENT) && (ORIG_SED)
Use Mod_Par
Use Mod_Prec 
Use Mod_Types
Use Mod_wd
Use Control, only : seddis
Use all_vars,only : CNSTNT,UNIFORM,SEDIMENT_PARAMETER_TYPE
# if defined (WAVE_CURRENT_INTERACTION)
USE MOD_WAVE_CURRENT_INTERACTION
# endif

# if defined (FLUID_MUD)
USE MOD_FLUID_MUD
# endif

implicit none 

!--------------------------------------------------
!Sediment Type                         
!
! sname        => sediment name (silt,clay,etc)     
! stype        => sediment type: 'cohesive'/'non-cohesive'
! Sd50         => sediment mean diameter (mm)
! Wset         => mean sediment settling velocity (mm/s)
! tau_ce       => critical shear stress for erosion (N/m^2)
! tau_cd       => critical shear stress for deposition 
! Srho         => sediment density (kg/m^3)
! Spor         => sediment porosity (dimensionless)
! erate        => surface erosion mass flux         [kg m^-2 s^-1] 
! conc         => sed concentration in water column [kg m^-3]
! cnew         => sed concentration during update   [kg m^-3]
! mass         => sediment mass in bed layers       [kg m^-2]
! frac         => sediment fraction in bed layers   [-]
! bflx         => bedload sediment flux (kg/m^2)  (+out of bed) 
! eflx         => suspended sediment erosive flux (+out of bed)
! dflx         => suspended sed depositional flux (+into bed)
! cdis         => concentration at river source 
! cflx         => store advective flux at open bndry 
! cobc         => user specd open bndry concentration 
! depm         => store deposited mass
! arraysize    => spatial size of sediment arrays
!--------------------------------------------------
type sed_type
  character(len=20) :: sname    
  character(len=20) :: sname2    
  character(len=20) :: stype 
  real(sp)          :: Sd50
  real(sp)          :: Wset  
  real(sp)          :: w0   ! settling velocity without floccuation
  real(sp)          :: tau_ce
  real(sp)          :: tau_cd
  real(sp)          :: Srho 
  real(sp)          :: Spor 
! Jianzhong Ge 03/05/2013
  real(sp)          :: Chin ! concentration for hindered settling 
  real(sp)          :: Wrdc ! reduction scale for hindered settling velocity
! Jianzhong Ge 03/05/2013
  real(sp)          :: erate 
  real(dp)          :: cmax 
  real(dp)          :: cmin 
  real(dp)          :: crms 
  real(sp), allocatable :: conc(:,:)
  real(sp), allocatable :: cnew(:,:)
  real(sp), allocatable :: mass(:,:)
  real(sp), allocatable :: frac(:,:)
  real(sp), allocatable :: bflx(:)        
  real(sp), allocatable :: eflx(:)        
  real(sp), allocatable :: dflx(:)        
  real(sp), allocatable :: cdis(:)
  real(sp), allocatable :: cflx(:,:)
  real(sp), allocatable :: cobc(:)
  real(sp), allocatable :: depm(:)   
! Jianzhong Ge 03/05/2013
  real(sp), allocatable :: wset0(:,:) ! vertical settling velocity
! Jianzhong Ge 03/05/2013
!-------------------Jianzhong Ge-------------------
  real(sp), allocatable :: t_cd(:)   ! =tau_cd
  real(sp), allocatable :: t_ce(:)   ! =tau_ce
  real(sp), allocatable :: rate(:)   ! =erate
!--------------------------------------------------
  integer           :: arraysize
end type
 
public

!--------------------------------------------------
!Global Model Parameters
!
! sedfile   : sediment input parameter control file
! nsed      : number of sediment classes
! nbed      : number of layers in sediment bed
! min_Srho  : minimum Sediment density
! inf_bed   : true if bed has infinite sediment supply
! bedload   : true if bedload is to be considered
! susload   : true if suspended load is to be considered
! DTsed     : sediment model time step (seconds)
! T_model   : model time (seconds)
! taub      : array to hold bottom shear stress
! rho0      : mean density parameter used to convert tau/rho to tau
! thck_cr   : critical thickness for initiating new surface layer (m)
! n_report  : iteration interval for statistics printing
! sed_start : start interval for sed model
! sed_its   : sediment model iteration counter
! sed_nudge : flag for activiating sediment nudging on obc
! sed_alpha : sediment nudging relaxation factor 
! sed_ramp  : number of iterations over which to ramp sed nudging
! sed_source: flag for activiating sediment point sources
!--------------------------------------------------
character(len=120), parameter :: sed_model_version = "cstms 1.0: non-cohesive only"
character(len=120) :: sedfile
integer  :: nsed
integer  :: nbed
integer  :: Nobc
integer  :: Nobc_gl
real(sp) :: min_Srho
logical  :: inf_bed
logical  :: sed_nudge
real(sp) :: sed_alpha
integer  :: sed_ramp  
logical  :: sed_dumpbed
logical  :: sed_dumpbot 
logical  :: sed_source
!Jianzhong Ge 03/05/2013
logical  :: vert_hindered
!Jianzhong Ge 03/05/2013
logical  :: morpho_model = .false.
real(sp) :: morpho_factor = 1.0
integer  :: morpho_incr = 1
integer  :: morpho_strt = 0
real(sp) :: rho_water = 1025.
logical  :: bedload
logical  :: oned_model
logical  :: susload
integer  :: n_report
integer  :: sed_start
logical  :: sed_hot_start
integer  :: sed_its
real(sp) :: DTsed
real(sp) :: T_model
real(sp) :: settle_cfl = 1.0
integer  :: settle_limiter = 2   
real(sp) :: qlim  = 1.0
real(sp), allocatable :: taub(:)
real(sp), parameter   :: rho0 = 1025.
real(sp), parameter   :: thck_cr = .005 
!--------------Jianzhong Ge------------------------
real(sp), parameter   :: mf=0.012,Vcons=1.e-6
real(sp), allocatable :: csed(:,:)
real(sp) :: deposition_d50
real(sp) :: sed_rhos
real(sp) :: sed_wset
!--------------------------------------------------
logical, parameter    :: debug_sed = .false.

!--------------------------------------------------
!Initial condition parameters
!  Read in sediment input file
!  Can be overwritten in init_sed.F or through 
!  restart
!
!  init_bed_porosity:  initial bed porosity [0,1]
!  init_bed_thickness: initial bed thickness [m] [> 0]
!  init_bed_fraction:  initial bed fraction [-] [> 0]
!--------------------------------------------------
 
real(sp) :: init_bed_porosity  = 1.0
real(sp), allocatable :: init_bed_thickness(:) 
real(sp), allocatable :: init_bed_fraction(:) 

!--------------------------------------------------
!Bedload 
!
!Meyer-Peter Muller Bedload Formulation Parameters
!
! Shield_Cr_MPM => Effective Critical Shields number
! Gamma_MPM     => MPM Power Law Coefficient
! k_MPM         => MPM Multiplier
!
!Generic Params for all bedload
!
! bedload_rate  => bedload rate coefficient
!
!--------------------------------------------------

real(sp) :: Shield_Cr_MPM  = 0.047  !  default: 0.047 
real(sp) :: Gamma_MPM      = 1.5    !  default: 1.5
real(sp) :: k_MPM          = 8.0    !  default: 8.0
real(sp) :: bedload_rate   = 0.1    !  default: 0.1
logical  :: bedload_smooth = .false.!  default: 0.1


!--------------------------------------------------
!Bed 
! data --> bed(horizontal,vertical,proptype)
! size --> bed(0:mt , nbed , n_bed_chars))
!
! with proptype=:
!   ithck  => layer thickness
!   iaged  => layer age 
!   iporo  => layer porosity  [volume of voids/total volume]
!
! global bed characteristics 
!
! n_bed_chars => number of bed characteristics 
!-------------------------------------------------

integer, parameter    :: n_bed_chars = 3
!character(len=80), dimension(n_bed_chars) :: bed_snames = [character(len=80) :: & 
!     "bed_thick","bed_age","bed_por"]
!character(len=80), dimension(n_bed_chars) :: bed_lnames = [character(len=80) :: & 
!     "bed thickness","bed age","bed porosity"]
!character(len=80), dimension(n_bed_chars) :: bed_units =  [character(len=80) :: & 
!     "m","days","-"]
character(len=80), dimension(n_bed_chars) :: bed_snames = (/"bed_thick","bed_age","bed_por"/)
character(len=80), dimension(n_bed_chars) :: bed_lnames = (/"bed thickness","bed age","bed porosity"/)
character(len=80), dimension(n_bed_chars) :: bed_units = (/"m","days","-"/)

integer, parameter    :: ithck = 1
integer, parameter    :: iaged = 2
integer, parameter    :: iporo = 3
real(sp), allocatable, target :: bed(:,:,:)

!--------------------------------------------------
!Bottom (Exposed Sediment Layer)
! data --> bed(0:mt,n_bot_vars)
!
! with proptype=:
!   isd50  => mean grain diameter 
!   isphi  => mean grain diameter in phi values
!   idens  => mean grain density
!   iwset  => mean settling velocity
!   itauc  => mean critical erosion stress [Pa, N/m^2]
!   iactv  => active layer thickness            
!   nthck  => new total thickness of sediment layer
!   lthck  => last thickness of sediment layer 
!   dthck  => accumulated delta of layer thickness [m]
!   tmass  => total mass in sediment layer [kg/m^2]
!   morph  => differential in thickness of sediment layer
!
! global bottom characteristics 
!   n_bot_vars => number of bottom characteristics
!
!-------------------------------------------------


!  BOTTOM properties indices:                                          !
!                                                                      !
!   MBOTP    Number of bottom properties (array dimension).            !
!   idBott   IO indices for bottom properties variables.               !
!   isd50    Median sediment grain diameter (m).                       !
!   isphi    mean grain diameter in phi values (-)                     !
!   idens    Median sediment grain density (kg/m3).                    !
!   iwsed    Mean settling velocity (m/s).                             !
!   itauc    Mean critical erosion stress (m2/s2).                     !
!   irlen    Sediment ripple length (m).                               !
!   irhgt    Sediment ripple height (m).                               !
!   ibwav    Bed wave excursion amplitude (m).                         !
!   izdef    Default bottom roughness (m).                             !
!   izapp    Apparent bottom roughness (m).                            !
!   izNik    Nikuradse bottom roughness (m).                           !
!   izbio    Biological bottom roughness (m).                          !
!   izbfm    Bed form bottom roughness (m).                            !
!   izbld    Bed load bottom roughness (m).                            !
!   izwbl    Bottom roughness used wave BBL (m).                       ! 
!   iactv    Active layer thickness for erosive potential (m).         !
!   ishgt    Sediment saltation height (m).                            !
!   idefx    Erosion flux                                              !
!   idnet    Erosion or deposition                                     !
!   idoff    Offset for calculation of dmix erodibility profile (m)    !
!   idslp    Slope for calculation of dmix or erodibility profile      !
!   idtim    Time scale for restoring erodibility profile (s)          !
!   idbmx    Bed biodifusivity maximum                                 !
!   idbmm    Bed biodifusivity minimum                                 !
!   idbzs    Bed biodifusivity zs                                      !
!   idbzm    Bed biodifusivity zm                                      !
!   idbzp    Bed biodifusivity phi                                     !
!   idprp    Cohesive behavior                                         !
!   bflag    Bed Flag (=0, erosion/deposition disabled)                !
!                                                                      !
!=======================================================================


integer, parameter    :: n_bot_vars  = 22
!character(len=80), dimension(n_bot_vars) :: bot_snames = [character(len=80) :: &
!      "bot_sd50","bot_dens","bot_wset","bot_tauc","bot_actv","bot_nthck",&
!      "bot_lthck","bot_dthck","bot_tmass","morph","rip_length","&
!       rip_height","wav_excur","def_rough","app_rough","Niku_rough","b&
!       io_rough","bfm_rough","bld_rough","wave_bbl","bot_phi","bedflag"]
!
!character(len=80), dimension(n_bot_vars) :: bot_lnames = [character(len=80) :: &
!      "bot_sd50","bot_dens","bot_wset","bot_tauc","bot_actv","sedlayer_depth_current_iteration",& 
!      "sedlayer_dapth_last_iteration","sedlayer_depth_net_change","bot_tmass","morph","rip_length","&
!       rip_height","wav_excur","def_rough","app_rough","Niku_rough","b&
!       io_rough","bfm_rough","bld_rough","wave_bbl","bot_phi","bedflag"]
!
!character(len=80), dimension(n_bot_vars) :: bot_units = [character(len=80) :: &  
!       "m","kg/m^3","m/s","N/m^2","days","m","m","m","kg","m","m","m","m","m"&
!      ,"m","m","m","m","m","m","-","-"] 
character(len=80), dimension(n_bot_vars) :: bot_snames = &
      (/"bot_sd50","bot_dens","bot_wset","bot_tauc","bot_actv","bot_nthck",&
      "bot_lthck","bot_dthck","bot_tmass","morph","rip_length","&
      &rip_height","wav_excur","def_rough","app_rough","Niku_rough","b&
      &io_rough","bfm_rough","bld_rough","wave_bbl","bot_phi","bedflag"/)
character(len=80), dimension(n_bot_vars) :: bot_lnames = &
      (/"bot_sd50","bot_dens","bot_wset","bot_tauc","bot_actv","sedlayer_depth_current_iteration",&
      "sedlayer_dapth_last_iteration","sedlayer_depth_net_change","bot_tmass","morph","rip_length","&
      &rip_height","wav_excur","def_rough","app_rough","Niku_rough","b&
      &io_rough","bfm_rough","bld_rough","wave_bbl","bot_phi","bedflag"/)
character(len=80), dimension(n_bot_vars) :: bot_units = (/"m","kg&
     &/m^3","m/s","N/m^2","days","m","m","m","kg","m","m","m","m","m"&
     &,"m","m","m","m","m","m","-","-"/)

integer, parameter    :: isd50  = 1  
integer, parameter    :: idens  = 2
integer, parameter    :: iwset  = 3
integer, parameter    :: itauc  = 4
integer, parameter    :: iactv  = 5 
integer, parameter    :: nthck  = 6 
integer, parameter    :: lthck  = 7 
integer, parameter    :: dthck  = 8 
integer, parameter    :: tmass  = 9 
integer, parameter    :: morph  = 10
integer, parameter    :: irlen  = 11
integer, parameter    :: irhgt  = 12
integer, parameter    :: ibwav  = 13
integer, parameter    :: izdef  = 14
integer, parameter    :: izapp  = 15
integer, parameter    :: izNik  = 16
integer, parameter    :: izbio  = 17
integer, parameter    :: izbfm  = 18
integer, parameter    :: izbld  = 19
integer, parameter    :: izwbl  = 20
integer, parameter    :: isphi  = 21
integer, parameter    :: bflag  = 22

real(sp), allocatable, target :: bottom(:,:)

!--------------------------------------------------
!Sediment Point Source Data 
!   sbc_tm:  time map for source data
!   seddis:  source data
!--------------------------------------------------

!--------------------------------------------------
!Sediment Array                  
!--------------------------------------------------

type(sed_type), allocatable, target :: sed(:)

!J. Ge for tracer advection
type(sed_type), allocatable, target :: sed0(:),sed2(:)
!J. Ge for tracer advection

contains


!==========================================================================
! Allocate Data, Initialize Concentration Fields and Bed Parameters       
!==========================================================================
  Subroutine Setup_Sed(fname,Nin,Nin_gl,N_sed,N_sed_Max,Sed_Names) 
  use control, only : msr,casename,input_dir
  use lims,    only : m
  Use Mod_Utils, only: pstop
  implicit none
  integer, intent(in)  :: Nin,Nin_gl   
  character(len=80)    :: fname
  integer, intent(out) :: N_sed
  integer, intent(in ) :: N_sed_max   
  character(len=20)    :: Sed_Names(N_sed_max)
  !-Local--------------------------
  logical             :: fexist
  integer             :: ised
  integer             :: i,k

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: setup_sed" 

  Nobc = Nin
  Nobc_gl = Nin_gl
!----------------------------------------------------------------
! check arguments and ensure necessary files exist
!----------------------------------------------------------------

  !initialize sediment model iteration counter
  sed_its = 0

  !ensure sediment parameter file exists
  if(fname == "none" .or. fname == "NONE" .or. fname == "None")then
    sedfile = "./"//trim(input_dir)//"/"//trim(casename)//'_sediment.inp'
  else
    sedfile = "./"//trim(input_dir)//"/"//trim(fname)
  endif
    
  inquire(file=trim(sedfile),exist=fexist)
  if(.not.fexist)then
    write(*,*)'sediment parameter file: ',trim(sedfile),' does not exist'
    write(*,*)'stopping'
    call pstop
  end if

!----------------------------------------------------------------
! read in sediment parameter file                    
!----------------------------------------------------------------
  call read_sed_params

! J. Ge for fluid mud layer
# if defined (FLUID_MUD)
!----------------------------------------------------------------
! read in fluid mud layer parameter file
!----------------------------------------------------------------
  call initialize_fluid_mud
# endif
! J. Ge for fluid mud layer


!----------------------------------------------------------------
! allocate data                                      
!----------------------------------------------------------------
  call alloc_sed_vars

!----------------------------------------------------------------
! initialize fields                                  
!   these will be overwritten if this is a restart case
!----------------------------------------------------------------
  call init_sed 
!----------------------------------------------------------------
! setup open boundary condition nudging values       
!----------------------------------------------------------------
  if(sed_nudge) call setup_sed_obc 

!----------------------------------------------------------------
! setup point source river forcing data 
!----------------------------------------------------------------
!  if(numqbc_gl .and. sed_source) call setup_sed_PTsource
!  if(sed_source) call setup_sed_PTsource

  !--------------------------------------------------
  !Initialize Bed_Mass properties
  !--------------------------------------------------
   Do k=1,Nbed
     Do i=1,m
       Do ised=1,Nsed
         sed(ised)%mass(i,k) = bed(i,k,ithck)*       &
                              (1.0-bed(i,k,iporo))*  &
                               sed(ised)%Srho*       &
                               sed(ised)%frac(i,k)
       end do
     end do
  end do

  !--------------------------------------------------
  !Calculate Bottom Thicknesses
  !--------------------------------------------------

  do i=1,m
    do k=1,nbed
      bottom(i,nthck) = bottom(i,nthck) + bed(i,k,ithck)
    end do
  end do


!----------------------------------------------------------------
! report sediment parameters and statistics to screen
!----------------------------------------------------------------
  call report_sed_setup

!----------------------------------------------------------------
! convert parameters to MKS 
!----------------------------------------------------------------
  call convert_sed_params 

!----------------------------------------------------------------
! update bottom properties                           
!----------------------------------------------------------------
  call update_bottom_properties

!----------------------------------------------------------------
! set arguments for sed names and N_sed for FVCOM to use 
! in setting up river forcing
!----------------------------------------------------------------
  N_sed = Nsed
  do i=1,Nsed
    sed_names(i) = sed(i)%sname
  end do
  if(dbg_set(dbg_sbr)) write(ipt,*) "End: setup_sed" 

End Subroutine Setup_Sed

!=======================================================================
!Convert Sed Params from nonstandard (mm, etc) to mks units
!=======================================================================
  Subroutine Convert_Sed_Params 
  Implicit None
  integer :: i

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Convert_Sed_Params" 

  !convert settling velocity to m/s from input mm/s
  do i=1,nsed
    sed(i)%Wset = sed(i)%Wset*.001
    sed(i)%w0   = sed(i)%w0*0.001
  end do

  !convert mean grain diameter from mm to m 
  do i=1,nsed
    sed(i)%Sd50 = sed(i)%Sd50*.001
  end do

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Convert_Sed_Params" 

  End Subroutine Convert_Sed_Params

!=======================================================================
!Read Sediment Parameters From Sediment Input File
!=======================================================================
  Subroutine Read_Sed_Params
  Use Input_Util
  Implicit None
  integer linenum,i,k1,iscan
  real(sp)           :: ftemp
  character(len=120) :: stemp

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Read_Sed_Params " 

  !check if we are running 1D model
  Call Get_Val(oned_model,sedfile,'SED_ONED',line=linenum,echo=.false.)

  !read in number of sediment classes 
  Call Get_Val(nsed,sedfile,'NSED',line=linenum,echo=.false.)  

  !read in start interval for sed model 
  Call Get_Val(sed_start,sedfile,'SED_START',line=linenum,echo=.false.)  

  !hot start option
  Call Get_Val(sed_hot_start,sedfile,'SED_HOT_START',line=linenum,echo=.false.)  

  !read in morphology parameters
  Call Get_Val(morpho_model,sedfile,'MORPHO_MODEL',line=linenum,echo=.false.)
  Call Get_Val(morpho_factor,sedfile,'MORPHO_FACTOR',line=linenum,echo=.false.)  
  if(morpho_model)then
    Call Get_Val(morpho_incr,sedfile,'MORPHO_INCR',line=linenum,echo=.false.)  
    Call Get_Val(morpho_strt,sedfile,'MORPHO_STRT',line=linenum,echo=.false.)  
  else
    morpho_incr   = 1
  endif
  
  !deactivate morpho model - not ready
!  morpho_model = .false.
!  morpho_factor = 1.0

  !read in interation interval for reporting 
  Call Get_Val(n_report,sedfile,'N_REPORT',line=linenum,echo=.false.)  

  !read in number of bed layers 
  Call Get_Val(nbed,sedfile,'NBED',line=linenum,echo=.false.)

  !read in logical for infinite sediment supply 
  Call Get_Val(inf_bed,sedfile,'INF_BED',line=linenum,echo=.false.)

  !read in initial bed porosity  
  Call Get_Val(init_bed_porosity,sedfile,'INIT_BED_POROSITY',line=linenum,echo=.false.)
  if(init_bed_porosity < 0. .or. init_bed_porosity > 1)then
    write(*,*)'error in init_bed_porosity in sed param file'
    write(*,*)'you entered: ',init_bed_porosity
    write(*,*)'must be in range [0,1]'
  endif

  !read in initial bed thickness
  allocate(init_bed_thickness(nbed))
  if(nbed==1)then
    Call Get_Val(ftemp,sedfile,'INIT_BED_THICKNESS',line=linenum,echo=.false.)
    init_bed_thickness(1) = ftemp
  else
    open(unit=999,file=trim(sedfile),form='formatted')
    iscan = scan_file(999,"INIT_BED_THICKNESS",fvec = init_bed_thickness,nsze = nbed)
    if(iscan /= 0)then
      write(*,*)'problem reading INIT_BED_THICKNESS from sediment param file'
      stop
    endif
    close(999)
  endif
  if(minval(init_bed_thickness) <= 0.)then
    write(*,*)'error in init_bed_thickness in sed param file'
    write(*,*)'you entered: ',init_bed_thickness
    write(*,*)'must be > 0'
  endif

  !read in initial bed fraction   
  allocate(init_bed_fraction(nsed))
  if(nsed==1)then
    Call Get_Val(ftemp,sedfile,'INIT_BED_FRACTION',line=linenum,echo=.false.)
    init_bed_fraction(1) = ftemp
  else
    open(unit=999,file=trim(sedfile),form='formatted')
    iscan = scan_file(999,"INIT_BED_FRACTION",fvec = init_bed_fraction,nsze = nsed)
    if(iscan /= 0)then
      write(*,*)'problem reading INIT_BED_FRACTION from sediment param file'
      stop
    endif
    close(999)
  endif
  if(minval(init_bed_fraction) < 0. .or. maxval(init_bed_fraction) > 1.)then
    write(*,*)'error in init_bed_fraction in sed param file'
    write(*,*)'also error in your chosen career'
    write(*,*)'you entered: ',init_bed_fraction  
    write(*,*)'must be >= 0 and <= 1'
  endif

  !read in logical for bedload calculation 
  Call Get_Val(bedload,sedfile,'BEDLOAD',line=linenum,echo=.false.)

  !shut down bedload if model is oned
  if(oned_model) bedload = .false.

  !read in logical for suspended load calculation 
  Call Get_Val(susload,sedfile,'SUSLOAD',line=linenum,echo=.false.)

  !read in minumum sediment density
  Call Get_Val(min_Srho,sedfile,'MIN_SRHO',line=linenum,echo=.false.)

  !read in nudging switch 
  Call Get_Val(sed_nudge,sedfile,'SED_NUDGE',line=linenum,echo=.false.)

  !read in params controlling output 
  Call Get_Val(sed_dumpbed,sedfile,'SED_DUMPBED',line=linenum,echo=.false.)
  Call Get_Val(sed_dumpbot,sedfile,'SED_DUMPBOT',line=linenum,echo=.false.)

  !read in nudging relaxation factor 
  if(sed_nudge) then
  Call Get_Val(sed_alpha,sedfile,'SED_ALPHA',line=linenum,echo=.false.)
  Call Get_Val(sed_ramp ,sedfile,'SED_RAMP ',line=linenum,echo=.false.)
  endif

  !settling CFL
  Call Get_Val(settle_cfl ,sedfile,'SETTLE_CFL',line=linenum,echo=.false.)
  if(settle_cfl <= 0 .or. settle_cfl > 1.)then
    write(*,*)'error in mod_sed, settle_cfl must be in range [>0,1]'
    write(*,*)'you entered: ',settle_cfl
    stop
  endif

  !limiter type
  Call Get_Val(settle_limiter,sedfile,'SETTLE_LIMITER',line=linenum,echo=.false.)
  if(settle_limiter < 1 .or. settle_limiter > 3)then
    write(*,*)'error in mod_sed, settle_limiter must be [1,2,3]'
    write(*,*)'you entered: ',settle_limiter
    stop
  endif
  if(settle_limiter==1) qlim = 0.
  if(settle_limiter==2) qlim = 1.
  if(settle_limiter==3) qlim = 2.
 
  !hindered effect on vertical diffusion
  Call Get_Val(VERT_HINDERED,sedfile,'VERT_HINDERED',line=linenum,echo=.false.)

  !read in point source switch 
  Call Get_Val(sed_source,sedfile,'SED_PTSOURCE',line=linenum,echo=.false.)

  !check values
  if(nsed < 1)then
    write(*,*)'sediment input file must have at least one sediment class'
    write(*,*)'currently nsed = ',nsed
    stop
  endif
  if(nbed < 1)then
    write(*,*)'sediment input file must have at least one bed layer'
    write(*,*)'currently nbed = ',nbed
    stop
  endif
  if(morpho_factor < 0. )then
    write(*,*)'sediment morpho_factor must be positive'
    write(*,*)'currently is = ',morpho_factor
    stop
  endif
  if(morpho_incr < 0 )then
    write(*,*)'sediment morpho_incr must be positive'
    write(*,*)'currently is = ',morpho_incr  
    stop
  endif

  if(morpho_strt < 0 )then
    write(*,*)'sediment morpho_strt must be non-negative'
    write(*,*)'currently is = ',morpho_strt  
    stop
  endif

  !allocate sediment data space
  Allocate(sed(nsed))

!J. Ge for tracer advection
  Allocate(sed0(nsed))
  Allocate(sed2(nsed))
!J. Ge for tracer advection


  !read in sediment parameters
  k1 = 1

  do i=1,nsed

    !read SED_NAME and mark position
    Call Get_Val(stemp,sedfile,'SED_NAME',line=linenum,echo=.false.,start=k1)
    sed(i)%sname = stemp
    sed(i)%sname2 = trim(sed(i)%sname)//'_bload'
    k1        = linenum+1

    !read type
    Call Get_Val(stemp,sedfile,'SED_TYPE',line=linenum,echo=.false.,start=k1)
    sed(i)%stype = stemp

    !read mean diameter
    Call Get_Val(ftemp,sedfile,'SED_SD50',line=linenum,echo=.false.,start=k1)
!Jianzhong GE
    sed(i)%Sd50 = ftemp
    deposition_d50 = ftemp*0.001
!------------------Jianzhong---------------
!    sed(i)%w0 = -4.*4.27*Vcons/(ftemp*1.0e-3)/1.22+ &
!                sqrt((4.*4.27/1.22*Vcons/(ftemp*1.0e-3))**2+ &
!                4./(3.*1.22)*1.65*9.8*(ftemp*1.0e-3)  )
!------------------------------------------
    !read sediment density 
    Call Get_Val(ftemp,sedfile,'SED_SRHO',line=linenum,echo=.false.,start=k1)
    sed(i)%Srho = ftemp 
!Jianzhong GE
    sed_rhos = ftemp

    !read sediment settling rate 
    Call Get_Val(ftemp,sedfile,'SED_WSET',line=linenum,echo=.false.,start=k1)
    sed(i)%Wset = ftemp
    sed(i)%w0   = ftemp
!Jianzhong GE
    sed_wset = ftemp*0.001

    !read sediment surface erosion rate 
    Call Get_Val(ftemp,sedfile,'SED_ERAT',line=linenum,echo=.false.,start=k1)
    sed(i)%erate = ftemp

    !read sediment critical erosive shear stress 
    Call Get_Val(ftemp,sedfile,'SED_TAUE',line=linenum,echo=.false.,start=k1)
    sed(i)%tau_ce = ftemp

    !read sediment critical depositional shear stress 
    Call Get_Val(ftemp,sedfile,'SED_TAUD',line=linenum,echo=.false.,start=k1)
    sed(i)%tau_cd = ftemp

    !read sediment porosity
    Call Get_Val(ftemp,sedfile,'SED_PORS',line=linenum,echo=.false.,start=k1)
    sed(i)%Spor = ftemp

!Jianzhong Ge 03/05/2013
    if(VERT_HINDERED)then
      !read concentration for hindered settling  
      if(sed(i)%stype=='cohesive')then
        Call Get_Val(ftemp,sedfile,'SED_CHIN',line=linenum,echo=.false.,start=k1)
        sed(i)%Chin = ftemp
      end if

      !read reduction scale for settling velocity
      if(sed(i)%stype=='cohesive')then
        Call Get_Val(ftemp,sedfile,'SED_WRDC',line=linenum,echo=.false.,start=k1)
        sed(i)%Wrdc = ftemp
      end if
    end if
!Jianzhong Ge 03/05/2013

  end do

  ! read in bedload function parameters
  Call Get_Val(Shield_Cr_MPM,sedfile,'MPM_CS',line=linenum,echo=.false.)
  Call Get_Val(    Gamma_MPM,sedfile,'MPM_GM',line=linenum,echo=.false.)
  Call Get_Val(        k_MPM,sedfile,'MPM_K ',line=linenum,echo=.false.)
  Call Get_Val(bedload_rate ,sedfile,'BEDLOAD_RATE ',line=linenum,echo=.false.)
  Call Get_Val(bedload_smooth ,sedfile,'BEDLOAD_SMOOTH',line=linenum,echo=.false.)
  
  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Read_Sed_Params " 
 
  End Subroutine Read_Sed_Params 
     
!=======================================================================
! Advance Sediment Model by time step DTSED 
!=======================================================================
  Subroutine Advance_Sed(DTin,Tin,taub_in)  
  Use Input_Util
  Use Scalar
  Use Control, only : ireport,iint,msr,par
  Use Lims,    only : m,mt,nt,kbm1,numqbc,kb,nprocs,myid
  Use All_Vars,only : BACKWARD_ADVECTION,BACKWARD_STEP
!  Use Mod_OBCS,only : iobcn
# if defined (MULTIPROCESSOR)
  Use Mod_Par
# endif

  Implicit None
  real(sp), intent(in ) :: DTin,Tin
  real(sp), intent(in ) :: taub_in(m)

  integer :: i,k,ised,l1,l2,ierr,d_cdis,d_cflx
  character(len=4) :: fnum
  real(sp) :: fact,ufact,temp

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Advance_Sed "

  !------------------------------------------------------
  !Check for Initialization
  !------------------------------------------------------
  if(iint < sed_start) return

  !------------------------------------------------------
  !Increment Sed Model Counter and Report Sed Model Init 
  !------------------------------------------------------
  sed_its = sed_its + 1
  if(msr .and.  sed_its == 1)then
    write(*,*)'========Sed Model Initiated======='
  endif

  !------------------------------------------------------
  !set bed age to current model time if first sed call 
  !------------------------------------------------------
  if(sed_its ==1)then
    bed(:,:,iaged) = T_model
  endif

  !------------------------------------------------------
  !Set up Model Forcing  -> time step and bottom shear 
  !Convert taux/tauy to node-based values
  !Convert Model tau/rho to tau [N/m^2]
  !------------------------------------------------------
  DTsed = DTin
  T_model = Tin
  taub(1:m) = taub_in(1:m)*rho_water

  !------------------------------------------------------
  !Calculate Active Layer Thickness --> bottom(:,iactv)
  !------------------------------------------------------
  call calc_active_layer

  !------------------------------------------------------
  !Vertical Sediment Dynamics --> Sediment 'Model'
  !------------------------------------------------------
  if(bedload)then
    !calculate bedload flux => sed(1:nsed)%bflx 
    !call calc_bedload_flux
    call calc_bedload_flux2

    !adjust surface bed properties
    call update_surface_bed('bedload')

    !add new surface layer if necessary
    call add_new_layer

    !update bed layer stratigraphy
    call update_stratigraphy

    !adjust bottom properties
    call update_bottom_properties
  endif

  if(susload)then
    !calculated depositional flux from settling
    call calc_deposition

    !calculate erosive flux
    call calc_erosion

! J. Ge for fluid mud layer

# if defined (FLUID_MUD)

!# if defined (WAVE_CURRENT_INTERACTION)
!    call fluid_mud_source_sink(DTsed,sed(1)%Srho,bottom(:,bflag),bed(:,1,ithck),bed(:,1,iporo),sed(1)%frac(:,1),bustrcwmax,bvstrcwmax,,sed(1)%t_ce,sed(1)%rate)
!# else
!    call fluid_mud_source_sink(DTsed,sed(1)%Srho,bottom(:,bflag),bed(:,1,ithck),bed(:,1,iporo),sed(1)%frac(:,1),wubot,wvbot,sed(1)%t_ce,sed(1)%rate)
!# endif
    call fluid_mud_source_sink(DTsed,sed(1)%Srho,bottom(:,bflag),bed(:,1,ithck),bed(:,1,iporo),sed(1)%frac(:,1),taub,sed(1)%t_ce,sed(1)%rate)

    call internal_step_fluid_mud

    do ised=1,nsed
      do i=1,m
# if defined (WET_DRY)
        if(iswetn(i)==1)then
# endif
          sed(ised)%conc(i,kbm1) = sed(ised)%conc(i,kbm1) + Entrainment(i)/(d(i)*dz(i,kbm1))
          sed(ised)%dflx(i) = dewatering(i)
          sed(ised)%eflx(i) = BedErosion(i) 
# if defined (WET_DRY)
        else
          sed(ised)%conc(i,kbm1) = 0.0_sp
          sed(ised)%dflx(i) = 0.0_sp
          sed(ised)%eflx(i) = 0.0_sp
        end if
# endif
      end do
    end do

# endif

    !adjust surface bed properties
    call update_surface_bed('susload')

    !add new surface layer if necessary
    call add_new_layer

    !update bed layer stratigraphy
    call update_stratigraphy

    !adjust bottom properties
    call update_bottom_properties
  endif
   
  !------------------------------------------------------
  !Horizontal Advection of Sediment Concentration 
  !------------------------------------------------------
  if(.not. oned_model)then
    d_cdis = 1
     if(numqbc > 0) d_cdis = numqbc
     d_cflx = 1
     if(Nobc > 0) d_cflx = Nobc
      
     do i=1,nsed
      
      !set discharge sediment concentrations
      if(numqbc > 0 .and. sed_source)then
        sed(i)%cdis(1:numqbc) = seddis(1:numqbc,i)
      endif
  
      !3D advection

!      call adv_scal(sed(i)%conc, &
!                      sed(i)%cnew, &
!                      d_cdis,      &
!                      sed(i)%cdis(1:numqbc), &
!                      d_cflx,      &
!                      sed(i)%cflx, &
!                      DTsed ,sed_source)
  

!J. Ge for tracer advection
      call adv_scal(sed(i)%conc,sed0(i)%conc,sed2(i)%conc, &
                      sed(i)%cnew, &
                      d_cdis,      &
                      sed(i)%cdis(1:numqbc), &
                      d_cflx,      &
                      sed(i)%cflx, &
                      DTsed ,sed_source)

!J. Ge for tracer advection


      call fct_sed(sed(i)%conc,sed(i)%cnew)
  
    end do
  else
   do i=1,nsed
     sed(i)%cnew = sed(i)%conc
   end do
  endif 


  !------------------------------------------------------
  !Vertical Diffusion of Sediment Concentrations
  !------------------------------------------------------

!Jianzhong GE 03/05/2013
  if(vert_hindered)call hindered_diffusion
!Jianzhong GE 03/05/2013

  do i=1,nsed
    call vdif_scal(sed(i)%cnew,DTsed) 
  end do

  !------------------------------------------------------
  !Exchange Concentration on Processor Boundaries
  !------------------------------------------------------
# if defined (MULTIPROCESSOR)
  if(par)then
    do i=1,nsed
      call aexchange(nc,myid,nprocs,sed(i)%cnew)
    end do
  endif
# endif

  !------------------------------------------------------
  ! Set Open Boundary Conditions (nudging to background)
  !------------------------------------------------------
  if(.not. oned_model)then
    if(sed_ramp == 0) sed_ramp = 1
    temp = min(float(sed_its/sed_ramp), 1.0)*sed_alpha
    if(.not.sed_nudge) temp = 0.0
    do i=1,nsed 
       call bcond_scal_OBC(sed(i)%conc, &
                           sed(i)%cnew, &
                           sed(i)%cflx, &
                           sed(i)%cobc, &
                           DTsed,       &
                           temp) 
    end do
  endif

  !------------------------------------------------------
  ! Point Source Boundary Conditions 
  !------------------------------------------------------
  if(sed_source .and. .not. oned_model)then
    do i=1,nsed
      call bcond_scal_PTsource(sed(i)%conc, &
                               sed(i)%cnew, &
                               sed(i)%cdis) 
    end do
  endif

  !------------------------------------------------------
  ! Update Concentration Variables to next time level
  !------------------------------------------------------
  do i=1,nsed
    sed(i)%conc = sed(i)%cnew
  end do

  !------------------------------------------------------
  ! Limit to Non-Negative (modify to Venkat Limiter) 
  !------------------------------------------------------
  do i=1,nsed
    sed(i)%conc = max(sed(i)%cnew,0.0)
  end do
  
  !-----------------------Jianzhong----------------------  
  csed = 0.0
  do i=1,nsed
    csed = csed + sed(i)%conc
  end do
  !------------------------------------------------------ 
  
  !------------------------------------------------------
  !Exchange Concentration on Processor Boundaries
  !------------------------------------------------------
# if defined (MULTIPROCESSOR)
  if(par)then
    do i=1,nsed
      call aexchange(nc,myid,nprocs,sed(i)%cnew)
    end do
    call aexchange(nc,myid,nprocs,taub)
  endif
# endif

!J. Ge for tracer advection
  IF(BACKWARD_ADVECTION==.TRUE.)THEN
    IF(BACKWARD_STEP==1)THEN
      do i=1,nsed
        SED0(i)%conc = sed(i)%cnew
      end do
    ELSEIF(BACKWARD_STEP==2)THEN
      do i=1,nsed
        SED2(i)%conc  = SED0(i)%conc 
        SED0(i)%conc  = sed(i)%cnew 
      end do
    ENDIF
  ENDIF
!J. Ge for tracer advection

  !------------------------------------------------------
  ! Update Depth Delta and Morphology Array
  ! this array is positive if material has accumulated,
  ! negative if eroded
  !------------------------------------------------------
  call update_thickness_delta

  !------------------------------------------------------
  ! Report 
  !------------------------------------------------------
  if(mod(sed_its,n_report)==0)call sed_report    

!  call layer_report(1)
  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Advance_Sed "

  End Subroutine Advance_Sed

!=======================================================================
! Allocate Sediment Variables                      
!=======================================================================
  Subroutine Alloc_Sed_Vars
  use lims,     only: nt,mt,kb,numqbc,kbm1
!  use mod_obcs, only: iobcn
  implicit none
  integer i,tmp1,tmp2,tmp3

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Alloc_Sed_Vars "

  !ensure arrays have nonzero dimension
  tmp1 = max(numqbc,1)
  tmp2 = max(Nobc+1,1)
  tmp3 = max(Nobc  ,1)

  !allocate suspended sediment arrays
  allocate(csed(0:mt,kb)) ; csed = 0.0

  do i=1,nsed
    sed(i)%arraysize = mt+1
    allocate(sed(i)%conc(0:mt,kb  ))      ; sed(i)%conc = 0.0
    allocate(sed(i)%cnew(0:mt,kb  ))      ; sed(i)%cnew = 0.0
    allocate(sed(i)%mass(0:mt,nbed))      ; sed(i)%mass = 0.0
    allocate(sed(i)%frac(0:mt,nbed))      ; sed(i)%frac = 0.0
    allocate(sed(i)%bflx(0:mt     ))      ; sed(i)%bflx = 0.0
    allocate(sed(i)%eflx(0:mt     ))      ; sed(i)%eflx = 0.0
    allocate(sed(i)%dflx(0:mt     ))      ; sed(i)%dflx = 0.0
    allocate(sed(i)%cdis(500      ))      ; sed(i)%cdis = 0.0
!JQI NOV2021
    allocate(sed(i)%cflx(tmp3,kbm1))      ; sed(i)%cflx = 0.0 ! KURT GLAESEMANN - remove +1 from dimension
!JQI NOV2021
    allocate(sed(i)%cobc(tmp3     ))      ; sed(i)%cobc = 0.0
    allocate(sed(i)%depm(0:mt     ))      ; sed(i)%depm = 0.0
! Jianzhong Ge 03/05/2013
    allocate(sed(i)%wset0(0:mt,kb  ))     ; sed(i)%wset0 = 0.0
    sed(i)%wset0 = sed_wset
! Jianzhong Ge 03/05/2013

!    if(SEDIMENT_PARAMETER_TYPE/=UNIFORM)then
    allocate(sed(i)%t_cd(0:mt     ))      ; sed(i)%t_cd = 0.0
    allocate(sed(i)%t_ce(0:mt     ))      ; sed(i)%t_ce = 0.0
    allocate(sed(i)%rate(0:mt     ))      ; sed(i)%rate = 0.0
!    end if

!J. Ge for tracer advection
!    if(backward_advection==.true.)then
      allocate(sed0(i)%conc(0:mt,kb  ))     ; sed0(i)%conc = 0.0
      allocate(sed2(i)%conc(0:mt,kb  ))     ; sed2(i)%conc = 0.0
!    endif
!J. Ge for tracer advection
  end do

  !allocate bottom shear stress array 
  allocate(taub(0:mt)) ; taub = 0.0

  !allocate bed data
  allocate(bed(0:mt , nbed , n_bed_chars)) ; bed = 0.0

  !allocate bottom data
  allocate(bottom(0:mt , n_bot_vars)) ; bottom = 0.0

  !allocate discharge data
!  write(*,*)'allocating discharage arrays: ',numqbc
!  pause
!  allocate(seddis(numqbc,10)) ; seddis = 0.0 

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Alloc_Sed_Vars "

  End Subroutine Alloc_Sed_Vars


!==========================================================================
! Report Sediment Setup To Screen 
!==========================================================================

  Subroutine Report_Sed_Setup
  use control, only: msr 
  implicit none
  integer :: i

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Report_Sed_Setup "

  ! echo sediment model parameters to screen
  if(msr)then

  write(*,*)
  write(*,*)'------------------ Sediment Model Setup --------------------'
  if(oned_model)then
    write(*,*)'!  model                 : oneD   '
    write(*,*)'!  1D model: bedload     : deactivated '
  endif
  write(*,*)'!  nbed                  :',nbed    
  write(*,*)'!  nsed                  :',nsed    
  write(*,*)'!  n_report              :',n_report
  if(sed_start > 0)then
    write(*,*)'!  sed_start             :',sed_start
  else
    write(*,*)'!  sed_start             : no delay'  
  endif
  if(sed_hot_start)then
    write(*,*)'!  sed_hot_start         :',sed_hot_start
  endif
    
  write(*,*)'!  min_Srho              :',min_Srho
  if(sed_nudge)then
    write(*,*)'!  open boundary nudging :  active'
    write(*,*)'!  nudging relax factor  :',sed_alpha 
    write(*,*)'!  # nudging ramp its    :',sed_ramp  
  else
    write(*,*)'!  open boundary nudging :  not active'
  endif
  if(sed_source)then
    write(*,*)'!  point source forcing  :  active'
  else
    write(*,*)'!  point source forcing  :  not active'
  endif
  if(morpho_model)then
    write(*,*)'!  morphodynamics        :  active'
    write(*,*)'!  morpho_incr           :',morpho_incr   
    write(*,*)'!  morpho_strt           :',morpho_strt   
  else
    write(*,*)'!  morphodynamics        :  not active'
  endif
  write(*,*)'!  morpho_factor         :',morpho_factor
  if(bedload)then
    write(*,*)'!  bedload dynamics      :  active'
    write(*,*)'!  Critical Shields      :',Shield_Cr_MPM
    write(*,*)'!  Bedload Exponent      :',Gamma_MPM     
    write(*,*)'!  Bedload Constant      :',k_MPM     
    write(*,*)'!  Bedload Rate Coeff.   :',bedload_rate
    if(bedload_smooth)then
      write(*,*)'!  Bedload Smoother      :  active'
    endif 
  else
    write(*,*)'!  bedload dynamics      :  not active'
  endif
  if(susload)then
    write(*,*)'!  susload dynamics      :  active'
  else
    write(*,*)'!  susload dynamics      :  not active'
  endif
  if(inf_bed)then
    write(*,*)'!  sediment supply       :  infinite'
  else
    write(*,*)'!  sediment supply       :  finite'
  endif
  write(*,*)'!  Settling CFL          :',settle_cfl
  write(*,*)'!'
  write(*,*)'!   class        type   Sd50    Wset   tau_ce   '//&
            'tau_cd  Erate   Spor    Srho'
  do i=1,nsed
    write(*,'(1X,A1,1X,A12,1X,A6,6(F7.4,1X),F8.2)')'!', &
              trim(sed(i)%sname),trim(sed(i)%stype(1:6)), &
              sed(i)%Sd50,sed(i)%Wset, &
              sed(i)%tau_ce,sed(i)%tau_cd,sed(i)%erate,sed(i)%Spor, &
              sed(i)%Srho  
  end do

  end if

  ! echo sediment model initial conditions to screen
  call sed_report  

  write(*,*)'------------------------------------------------------------'
  write(*,*)

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Report_Sed_Setup "

  End Subroutine Report_Sed_Setup

!==========================================================================
! Monitor Sediment Model: Compute and Write Statistics to Screen
!==========================================================================

  Subroutine Sed_Report  
  use control, only: msr ,par, mpi_fvcom_group
  use lims   , only: m,mt,kbm1, msrid
  implicit none
  real(dp) :: tau_min,tau_max,tmp1(3),tmp2(3),dthck_max
  integer :: i,dim1,dim2,ierr
#  if defined (FLUID_MUD)
  real(dp) :: temp1(5),temp2(5),mud_max,ua_max,va_max,ua_min,va_min
  REAL(DP), DIMENSION(3) :: SBUF,RBUF1,RBUF2,RBUF3
  real(dp) :: ESTOT,NSTOT
#  endif

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Sed_Report "

  !set up limits
  dim1 = mt 
  dim2 = kbm1
  
  !calculate max/min/rms concentrations
  do i=1,nsed
    sed(i)%cmax = maxval(sed(i)%conc(1:dim1,1:dim2))
    sed(i)%cmin = minval(sed(i)%conc(1:dim1,1:dim2))
    sed(i)%crms = sum(dble(abs(sed(i)%conc(1:dim1,1:dim2))))
    sed(i)%crms = sed(i)%crms/(dim1*dim2)
  end do

  dthck_max = maxval(abs(bottom(1:m,dthck)))
  tau_max   = maxval(taub(1:m))
  tau_min   = minval(taub(1:m))

#  if defined (FLUID_MUD)
  ESTOT = DBLE(NGL)
  NSTOT = DBLE(MGL)
  !calculate/max/min/rms mud layer
  mud_max = maxval(abs(D_FMLcm))
  ua_max = maxval(ua_fml)
  ua_min = minval(ua_fml)
  va_max = maxval(va_fml)
  va_min = minval(va_fml)

  SBUF(1)  = SUM(DBLE(UA_fml(1:n)))
  SBUF(2)  = SUM(DBLE(VA_fml(1:n)))
  SBUF(3)  = SUM(DBLE(D_FML(1:m)))

  RBUF1 = SBUF

#  endif
 
# if defined (MULTIPROCESSOR)
  IF(PAR)THEN
  do i=1,nsed
    tmp1(1:3) = (/sed(i)%cmax,sed(i)%cmin,sed(i)%crms/)
    CALL MPI_REDUCE(tmp1,tmp2,3,MPI_DP,MPI_MAX,MSRID-1,MPI_FVCOM_GROUP,IERR)
    if(msr)then
      sed(i)%cmax = tmp2(1)
      sed(i)%cmin = tmp2(2)
      sed(i)%crms = tmp2(3)
    endif
  end do

  tmp1(1:3) = (/dthck_max,tau_min,tau_max/)
  CALL MPI_REDUCE(tmp1(1),tmp2(1),1,MPI_DP,MPI_MAX,MSRID-1,MPI_FVCOM_GROUP,IERR)
  CALL MPI_REDUCE(tmp1(2),tmp2(2),1,MPI_DP,MPI_MIN,MSRID-1,MPI_FVCOM_GROUP,IERR)
  CALL MPI_REDUCE(tmp1(3),tmp2(3),1,MPI_DP,MPI_MAX,MSRID-1,MPI_FVCOM_GROUP,IERR)
  dthck_max = tmp2(1) 
  tau_min   = tmp2(2) 
  tau_max   = tmp2(3) 

#  if defined (FLUID_MUD)
  temp1(1:5) = (/mud_max,ua_max,ua_min,va_max,va_min/)
  CALL MPI_REDUCE(temp1(1),temp2(1),1,MPI_DP,MPI_MAX,MSRID-1,MPI_FVCOM_GROUP,IERR)
  CALL MPI_REDUCE(temp1(2),temp2(2),1,MPI_DP,MPI_MAX,MSRID-1,MPI_FVCOM_GROUP,IERR)
  CALL MPI_REDUCE(temp1(3),temp2(3),1,MPI_DP,MPI_MIN,MSRID-1,MPI_FVCOM_GROUP,IERR)
  CALL MPI_REDUCE(temp1(4),temp2(4),1,MPI_DP,MPI_MAX,MSRID-1,MPI_FVCOM_GROUP,IERR)
  CALL MPI_REDUCE(temp1(5),temp2(5),1,MPI_DP,MPI_MIN,MSRID-1,MPI_FVCOM_GROUP,IERR)
  mud_max = temp2(1)
  ua_max  = temp2(2)
  ua_min  = temp2(3)
  va_max  = temp2(4)
  va_min  = temp2(5)
  CALL MPI_REDUCE(SBUF,RBUF1,3,MPI_DP,MPI_SUM,MSRID-1,MPI_FVCOM_GROUP,IERR)
  IF(ISNAN(RBUF1(1)).or.ISNAN(RBUF1(2)).or.ISNAN(RBUF1(3)))THEN
     CALL FATAL_ERROR("fluid mud caclulation fails with unstable iteration","Please check the model")
   END IF
#  endif

  ENDIF
# endif

  ! write statistics to screen
  if(.not.msr)return
  write(*,*  )'====================Sediment Model Stats======================'
  write(*,*  )'!  quantity              :     avg(g/L)      max(g/L)      '
  do i=1,nsed
    write(*,100)'!',trim(sed(i)%sname),'    :',sed(i)%crms,sed(i)%cmax
  end do
  write(*,*  )'!  max sed thick change  :     ',dthck_max
  write(*,*  )'!  max bottom stress     :     ',tau_max
  write(*,*  )'!  min bottom stress     :     ',tau_min
#  if defined (FLUID_MUD)
  write(*,*  )'!  max mud thickness     :     ',mud_max
  write(*,*  )'!  max u-speed of mud    :     ',ua_max
  write(*,*  )'!  min u-speed of mud    :     ',ua_min
  write(*,*  )'!  max v-speed of mud    :     ',va_max
  write(*,*  )'!  min v-speed of mud    :     ',va_min
  write(*,*  )'!  avg u-speed of mud    :     ',RBUF1(1)/ESTOT
  write(*,*  )'!  avg v-speed of mud    :     ',RBUF1(2)/ESTOT
  write(*,*  )'!  avg thickness of mud  :     ',RBUF1(3)/NSTOT
#  endif  

  100 format(1x,a1,a20,a5,3f12.6)
  101 format(1x,a26,f12.6)

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Sed_Report "

  End Subroutine Sed_Report  

!==========================================================================
! Calculate Bed Load Transport using Meyer-Peter Muller Formulation   
!   Modify for influence of bathy slope
!   See eqns 24-27+35 + Warner et al, Comp. Geos.  2008 
!==========================================================================

  Real(sp) Function Bedload_MPM(rho_sed,D_50,tau_b) 
  use control, only: grav
  implicit none
  real(sp), intent(in) :: rho_sed
  real(sp), intent(in) :: D_50 
  real(sp), intent(in) :: tau_b  
  !--------------------------------
  real(sp) :: tau_nd,bedld_nd,R

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Bedload_MPM "

  R = (rho_sed/1025.-1)

  !nondimensional bottom stress 
  tau_nd = tau_b/(R*grav*D_50) 

  !compute bedload flux (PHI, eq. 25)
  bedld_nd = k_MPM*(MAX((tau_nd - Shield_Cr_MPM),0.0)**Gamma_MPM)

  !compute bedload (kg/(m-s)) (q_bl eq. 24)
  bedload_MPM = sqrt(R*grav*D_50)*D_50*rho_sed*bedld_nd

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Bedload_MPM "

  End Function Bedload_MPM

!==========================================================================
! Given Bed Load Transports, Calculate the flux using divergence 
!   - calculate bedload transport in elements
!   - modify transport using bed slope
!   - exchange at interprocessor boundaries
!   - calculate divergence to give flux out of / into bed
!==========================================================================
  Subroutine  Calc_Bedload_Flux() 

  !use all_vars, only : niec,ntrg,ncv,nt,dltxe,dltye,nprocs,n,msr,par,wubot,wvbot,myid,art1,m
  use all_vars
  use bcs, only : iobcn,i_obc_n
  use mod_obcs, only: next_obc
# if defined (MULTIPROCESSOR)
  Use Mod_Par
# endif
  integer  :: i1,ia,ib,ns,i,n1,n2,n3
  real(sp) :: qbl,ootaub_e,taub_e,flux
  real(sp) :: dbdx,dbdy,beta,slopex,slopey,sed_slope,max_slope,bld_dr
  real(sp), parameter :: sed_angle = 33.0  !sed "friction angle" in degrees
  real(sp), allocatable, target :: btx(:) 
  real(sp), allocatable, target :: bty(:) 
  real(sp), parameter :: taub_min = tiny(1.0_sp)

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Calc_Bedload_Flux "

  !initialize working bedload transport arrays
  allocate(btx(0:nt))
  allocate(bty(0:nt))

  !main loop over sed types (cohesive?)
  do ns = 1,nsed
    sed(ns)%bflx = 0.0
    btx = 0.0
    bty = 0.0

    !loop over elements, calculate bedload transport in x (btx) and  y (bty) at each
    ! => btx (qblx, eq 28) [after loop is kg/m]
    ! => bty (qbly, eq 29) [after loop is kg/m]
    do i=1,n
  
      !magnitude (units stress/density) of bottom shear stress at elements
      taub_e = sqrt(wubot(i)**2 + wvbot(i)**2)
    
      !calculate transport using Meyer-Peter Mueller formulation
      !qbl has units of kg-m^2/s
      qbl = bedload_mpm(sed(ns)%Srho,sed(ns)%Sd50,taub_e)
      if(taub_e > taub_min)then
        ootaub_e = 1./taub_e
      else
        ootaub_e = 0.0
      endif


      !calculate the bedload transport (kg/m)
      !note the negative sign: bottom stress in FVCOM (wubot/wvbot) is stress
      !on the fluid.  So a + current will produce a - bed stress
      !we want this transport to be in the direction of the current, thus the 
      !change of sign
      btx(i) = -DTsed*qbl*wubot(i)*ootaub_e  
      bty(i) = -DTsed*qbl*wvbot(i)*ootaub_e 
  
    end do
  
    !add bed slope contribution to our horizontal fluxes
    !Eq 35, Warner etal, Comp & Geo, 2008 => qbl_slope = slopex,slopey 
    !minimize the maximum downslope, but do not truncate upslope
    !upslope is limited by Lesser formula
    !range of slopefactor (slopex)  = [.55,5.64] for ..
    !..  upslope(dhdx<0)=-infty -> downslope(dhdx>0) = infty
    !dbdx > 0 is a downslope (see bld_dr to enforce this)
    sed_slope = tan(sed_angle*3.14159/180.)
    max_slope = 0.8*sed_slope
    do i=1,n
      n1 = nv(i,1) ; n2 = nv(i,2) ; n3 = nv(i,3)

      bld_dr = sign(1.0_sp,btx(i))
      dbdx   = min( (awx(i,1)*h(n1)+awx(i,2)*h(n2)+awx(i,3)*h(n3))*bld_dr ,max_slope)
      beta   = atan(dbdx)
      slopex = sed_slope/( (sed_slope-dbdx)*cos(beta)) 
      btx(i) = btx(i)*slopex

      bld_dr = sign(1.0_sp,bty(i))
      dbdy   = min( (awy(i,1)*h(n1)+awy(i,2)*h(n2)+awy(i,3)*h(n3))*bld_dr ,max_slope)
      beta   = atan(dbdy)
      slopey = sed_slope/( (sed_slope-dbdy)*cos(beta)) 
      bty(i) = bty(i)*slopey
      
    end do
  
    ! multiply by morphological time scale (to speed up time)
!    if(morpho_model)then
       btx = btx*morpho_factor
       bty = bty*morpho_factor
!    endif
  
    ! multiply by bedload transport coefficient (fudge factor) 
     btx = btx*bedload_rate 
     bty = bty*bedload_rate 
  
    ! limit transport to available bed mass
    ! ROMS does this, but doesn't make sense, mass is linked to divergence of  
    ! transport, not horizontal fluxes , check with J. Warner
     
  
    !exchange btx/bty across interprocessor boundaries 
# if defined (MULTIPROCESSOR)
    if(par)then
      call aexchange(ec,myid,nprocs,btx,bty)
    endif
# endif
  
    !calculate the divergence to give flux (kg) on the nodes
    do i=1,ncv
       i1   = ntrg(i)
       ia   = niec(i,1)
       ib   = niec(i,2)
       flux = (-btx(i1)*dltye(i) + bty(i1)*dltxe(i))
       sed(ns)%bflx(ia) = sed(ns)%bflx(ia)-flux
       sed(ns)%bflx(ib) = sed(ns)%bflx(ib)+flux
     end do

     !update the final bedload flux (kg/m^2) on the nodes
     !this flux is + out of the bed
     !if flow dir is positive and bed_stress is increasing in x
     !this should produce a net positive (out of bed) flux
     do i=1,m
       sed(ns)%bflx(i) = sed(ns)%bflx(i)/art1(i)
     end do

     !set the flux on the open boundary
     do i=1,iobcn
       sed(ns)%bflx( i_obc_n(i) ) = sed(ns)%bflx(next_obc(i))
     end do

     !shutdown bedload if node is locally an inactive node
     do i=1,m
       sed(ns)%bflx(i) = sed(ns)%bflx(i)*bottom(i,bflag)
     end do

     !if bathy is rough - choose smooth_bedload, averaging filter
     if(bedload_smooth)then
       call n2e2d(sed(ns)%bflx,btx)
       call e2n2d(btx,sed(ns)%bflx)
     endif
       
  end do !loop over sediment types

  deallocate(btx)
  deallocate(bty)

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Calc_Bedload_Flux "  

  End Subroutine Calc_Bedload_Flux

!==========================================================================
! Given Bed Load Transports, Calculate the flux using divergence 
!   - calculate bedload transport in elements
!   - modify transport using bed slope
!   - exchange at interprocessor boundaries
!   - calculate flux on elements, upwinding using sign of tau dot n 
!   - convert bedload flux to nodes for morphology
!==========================================================================
  Subroutine  Calc_Bedload_Flux2() 

  !use all_vars, only : niec,ntrg,ncv,nt,dltxe,dltye,nprocs,n,msr,par,wubot,wvbot,myid,art1,m
  use all_vars
  use bcs, only : iobcn,i_obc_n
  use mod_obcs, only: next_obc
# if defined (MULTIPROCESSOR)
  Use Mod_Par
# endif
  integer  :: i1,ia,ib,ns,i,n1,n2,n3,j1,j2,cnt,jj,cc
  real(sp) :: qbl,ootaub_e,taub_e,flux
  real(sp) :: dbdx,dbdy,beta,slopex,slopey,sed_slope,max_slope,bld_dr
  real(sp) :: tau_x,tau_y,flowdir,bldx,bldy
  real(sp), parameter :: one = 1.0_sp
  real(sp), parameter :: eps = .5    
  real(sp), parameter :: sed_angle = 33.0  !sed "friction angle" in degrees
  real(sp), allocatable, target :: btx(:) 
  real(sp), allocatable, target :: bty(:) 
  real(sp), allocatable, target :: bflxe(:) 
  real(sp), parameter :: taub_min = tiny(1.0_sp)

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Calc_Bedload_Flux2 "

  !initialize working bedload transport arrays
  allocate(btx(0:nt))
  allocate(bty(0:nt))
  allocate(bflxe(0:nt))

  !main loop over sed types (cohesive?)
  do ns = 1,nsed
    sed(ns)%bflx = 0.0
    bflxe = 0.0
    btx   = 0.0
    bty   = 0.0

    !loop over elements, calculate bedload transport in x (btx) and  y (bty) at each
    ! => btx (qblx, eq 28) [after loop is kg/m]
    ! => bty (qbly, eq 29) [after loop is kg/m]
    do i=1,n
  
      !magnitude (units stress/density) of bottom shear stress at elements
      taub_e = sqrt(wubot(i)**2 + wvbot(i)**2)
    
      !calculate transport using Meyer-Peter Mueller formulation
      qbl = bedload_mpm(sed(ns)%Srho,sed(ns)%Sd50,taub_e)
      if(taub_e > taub_min)then
        ootaub_e = 1./taub_e
      else
        ootaub_e = 0.0
      endif


      !calculate the bedload transport (kg/m)
      !note the negative sign: bottom stress in FVCOM (wubot/wvbot) is stress
      !on the fluid.  So a + current will produce a - bed stress
      !we want this transport to be in the direction of the current, thus the 
      !change of sign
      btx(i) = -DTsed*qbl*wubot(i)*ootaub_e  
      bty(i) = -DTsed*qbl*wvbot(i)*ootaub_e 
  
    end do
  
    !add bed slope contribution to our horizontal fluxes
    !Eq 35, Warner etal, Comp & Geo, 2008 => qbl_slope = slopex,slopey 
    !minimize the maximum downslope, but do not truncate upslope
    !upslope is limited by Lesser formula
    !range of slopefactor (slopex)  = [.55,5.64] for ..
    !..  upslope(dhdx<0)=-infty -> downslope(dhdx>0) = infty
    !dbdx > 0 is a downslope (see bld_dr to enforce this)
    sed_slope = tan(sed_angle*3.14159/180.)
    max_slope = 0.8*sed_slope
    do i=1,n
      n1 = nv(i,1) ; n2 = nv(i,2) ; n3 = nv(i,3)

      bld_dr = sign(1.0_sp,btx(i))
      dbdx   = min( (awx(i,1)*h(n1)+awx(i,2)*h(n2)+awx(i,3)*h(n3))*bld_dr ,max_slope)
      beta   = atan(dbdx)
      slopex = sed_slope/( (sed_slope-dbdx)*cos(beta)) 
      btx(i) = btx(i)*slopex

      bld_dr = sign(1.0_sp,bty(i))
      dbdy   = min( (awy(i,1)*h(n1)+awy(i,2)*h(n2)+awy(i,3)*h(n3))*bld_dr ,max_slope)
      beta   = atan(dbdy)
      slopey = sed_slope/( (sed_slope-dbdy)*cos(beta)) 
      bty(i) = bty(i)*slopey
      
    end do
  
    ! multiply by morphological time scale (to speed up time)
!    if(morpho_model)then
       btx = btx*morpho_factor
       bty = bty*morpho_factor
!    endif
  
    ! multiply by bedload transport coefficient (fudge factor) 
     btx = btx*bedload_rate 
     bty = bty*bedload_rate 
  
    ! limit transport to available bed mass
    ! ROMS does this, but doesn't make sense, mass is linked to divergence of  
    ! transport, not horizontal fluxes , check with J. Warner
     
  
    !exchange btx/bty across interprocessor boundaries 
# if defined (MULTIPROCESSOR)
     if(par)then
       call aexchange(ec,myid,nprocs,btx,bty)
     endif
# endif

     !set dummy values for btx/bty to check divergence calc
     !do i=1,n
     !  btx(i) = 0.0 !-xc(i)
     !  bty(i) = -yc(i) 
     !end do
  
     !calculate the flux on the elements using first order upwind formulation 
     do i=1,ne
       ia=iec(i,1)
       ib=iec(i,2)
       j1=ienode(i,1)
       j2=ienode(i,2)
       !bed stress in x and y averaged to an edge , used for upwind selection
       tau_x = 0.5*(wubot(ia)+wubot(ib))  
       tau_y = 0.5*(wvbot(ia)+wvbot(ib))  
       !dot the bed stress vector with length-weighted edge normal 
       !note bed stress opposes flow dir so we need neg sign
       !the vector (-dltyc,dltxc) points from ia to ib
       flowdir  = -(-tau_x*dltyc(i) + tau_y*dltxc(i)) ! -(tau dot dL)
       !upwind select the x-dir and y-dir bedload transport 
       bldx  = ((one-sign(one,flowdir))*btx(ib)+(one+sign(one,flowdir))*btx(ia))*0.5_sp 
       bldy  = ((one-sign(one,flowdir))*bty(ib)+(one+sign(one,flowdir))*bty(ia))*0.5_sp 
       !accum contributions to bedload transport divergence calculations at elements ia/ib
       !on either side of edge i
       bflxe(ia) = bflxe(ia) - bldx*dltyc(i) + bldy*dltxc(i)
       bflxe(ib) = bflxe(ib) + bldx*dltyc(i) - bldy*dltxc(i)
     end do

     !divide by element area to complete divergence => bedload in kg/m^2
     !if there was a divergence of bedload transport, this should be a positive number
     do i=1,n
       bflxe(i) = bflxe(i)/art(i)
     end do
      
     !exchange bflxe across interprocessor boundaries
#    if defined (MULTIPROCESSOR)
     if(par)then
       call aexchange(ec,myid,nprocs,bflxe)
     endif
#    endif

     !interpolate element-based blfx to the nodes
     call e2n2d(bflxe,sed(ns)%bflx)

     !smooth using explicit smoother
     do i=1,m
       do jj=1,ntsn(i)
         cc = nbsn(i,jj)
         sed(ns)%bflx(i) = sed(ns)%bflx(i) + eps*sed(ns)%bflx(cc)
       end do
       sed(ns)%bflx(i) = sed(ns)%bflx(i)/(1. + eps*ntsn(i))
     end do

     !set the flux on the open boundary
     do i=1,iobcn
       sed(ns)%bflx( i_obc_n(i) ) = sed(ns)%bflx(next_obc(i))
     end do

     !shutdown bedload if node is locally an inactive node
     do i=1,m
       sed(ns)%bflx(i) = sed(ns)%bflx(i)*bottom(i,bflag)
     end do

  end do !loop over sediment types

  !cleanup 
  deallocate(btx)
  deallocate(bty)
  deallocate(bflxe)
  
  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Calc_Bedload_Flux2 "

  End Subroutine Calc_Bedload_Flux2



!==========================================================================
! Calculate Erosion Rate (kg/m^2) 
! Calculate erosion rate using the method of Ariathurai and Arulanandan (1978)
!    -Erosion rates of cohesive soils. J. Hydraulics Division, ASCE, 104(2),279-282.-
!    surface_mass_flux[i] = erate[i]*(1-porosity)*bfrac[i]*(tau_w/tau_crit[i] - 1)
!      where
!      i               = index of sediment type
!      bfrac[i]        = fraction of bed composed by sediment i
!      erate[i]        = bed erodibility constant
!      porosity        = volume of voids/total volume in the top layer
!      tau_w           = shear stress on the bed
!      tau_crit[i]     = critical shear stress of 
!==========================================================================

  Subroutine Calc_Erosion
  use all_vars
  implicit none
  integer  :: ised ,i
  real(sp) :: bed_por,bed_frac,tau_w,tau_ce,erate,dep,active,erosion_orig
  real(sp) :: dx,dy 
  integer  :: cnt,jj,cc
  real(sp),dimension(mt) :: unode,vnode

  ! ------ new: Karsten Lettmann ---
  real(sp) :: dz_bed, eflux_avail
  ! --------------------------------

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Calc_Erosion "
  ! Compute bottom velocity at node
   if(.not.oned_model)then
   CALL E2N2D(u(:,kbm1),unode)
   CALL E2N2D(v(:,kbm1),vnode)
   endif

  !Compute Erosive Flux (kg/m^2), limited by available material in active layer 
  do ised=1,nsed
    do i=1,m
# if defined (WET_DRY)
      if(iswetn(i)/=1)cycle  ! no calculation when dry node
# endif      
      if(SEDIMENT_PARAMETER_TYPE/=UNIFORM)then
        erate  = sed(ised)%rate(i)
        tau_ce = sed(ised)%t_ce(i) 
      else
        erate  = sed(ised)%erate
        tau_ce = sed(ised)%tau_ce 
      end if 
      !calc taub from upwind or downwind elements - beta
# if defined (WAVE_CURRENT_INTERACTION)
      if(MB_BBL_USE.or.SG_BBL_USE.or.SSW_BBL_USE)then
         taub(i) = 1025.0 * taucwmax(i)
      else
# endif
        taub(i) = 0.
        cnt = 0.
        do jj=1,ntve(i)
          cc = nbve(i,jj)
          !if(xc(cc) < vx(i))then !cell is upstream, only for eastward flow
          !if(xc(cc) > vx(i))then !cell is downtream, only for eastward flow
          !dx=xc(cc)-vx(i);dy=yc(cc)-vy(i);
          !if((dx*unode(i)-dy*vnode(i))>0.) then ! arbitrary direction
            cnt = cnt + 1
            taub(i) = taub(i) + sqrt(wubot(cc)**2 + wvbot(cc)**2) 
          !endif
        end do
        taub(i) = 1025.*taub(i)/cnt
# if defined (WAVE_CURRENT_INTERACTION)
      end if
# endif
     
! J. Ge for fluid mud layer
# if !defined (FLUID_MUD)
      bed_por  = bed(i,1,iporo)
      bed_frac = sed(ised)%frac(i,1)
      tau_w    = taub(i)
      dep      = sed(ised)%dflx(i)
      active   = (1.0-bed_por)*Bed_Frac*(tau_w/tau_ce-1.0)*bottom(i,iactv) + dep
      !erosion_orig  = MAX(DTsed*Erate*(1.0-bed_por)*Bed_Frac*(tau_w/tau_ce-1.0),0.0)
      !sed(ised)%eflx(i) = erosion_orig 
      sed(ised)%eflx(i) = DTsed*erosion(erate,bed_por,bed_frac,tau_w,tau_ce) !min(erosion,active) !geoff
      !sed(ised)%eflx(i) = DTsed*min(erosion(erate,bed_por,bed_frac,tau_w,tau_ce),active) !min(erosion,active) !geoff

      !shutdown erosion if node is inactive 
      sed(ised)%eflx(i) = sed(ised)%eflx(i)*bottom(i,bflag)

      ! ========= new: Karsten Lettmann, Okt. 2012 ==========
      ! limit the eflux by the thickness of the top layer_report
      ! how much is available for erosion
      dz_bed = bed(i,1,ithck)
      eflux_avail = bed_frac * sed(ised)%srho*(1.0-bed_por) * dz_bed + dep

      ! take only 99% of possible material
      sed(ised)%eflx(i) = min(sed(ised)%eflx(i) , eflux_avail)
      ! ===============  end new ================================

!      if(i==1)print*,'calc_erosion:',dep,bed(i,1,ithck),eflux_avail,sed(ised)%eflx(i)

! J. Ge for fluid mud layer
# endif

    end do
  end do

! J. Ge for fluid mud layer

# if !defined (FLUID_MUD)
  !Update Concentration in Bottom of Water Column
  do ised=1,nsed
    do i=1,m
# if defined (WET_DRY)
      if(iswetn(i)==1)then
# endif
        sed(ised)%conc(i,kbm1) = sed(ised)%conc(i,kbm1) + sed(ised)%eflx(i)/(d(i)*dz(i,kbm1))
# if defined (WET_DRY)
      else
        sed(ised)%conc(i,kbm1) = 0.0_sp
      end if
# endif
    end do
  end do
# endif
  
  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Calc_Erosion "

  End Subroutine Calc_Erosion

!==========================================================================
! Calculate Erosion (kgm^-2s^-1) using user-defined formula 
!==========================================================================
  Real(sp) Function Erosion(erate,bed_por,bed_frac,tau_w,tau_ce)
  implicit none
  real(sp), intent(in) ::  erate 
  real(sp), intent(in) ::  bed_por 
  real(sp), intent(in) ::  bed_frac 
  real(sp), intent(in) ::  tau_w 
  real(sp), intent(in) ::  tau_ce 
!Jianzhong
  real(sp),parameter :: a1=-0.144,a2=0.904,a3=-0.823,a4=0.204  
  real(sp)::ratio
!JIanzhong

  !standard CSTM formulation 
!  Erosion = MAX(Erate*(1.0-bed_por)*Bed_Frac*(tau_w/tau_ce-1.0),0.0)

!Jianzhong
  !Prooijen-Winterwerp formulation
   ratio=tau_w/tau_ce
   if(ratio<0.52)then
     Erosion=0.0
   elseif(ratio>1.7)then
     Erosion = Erate*(1.0-bed_por)*Bed_Frac*(tau_w/tau_ce-1.0)    
   else
     Erosion = Erate*(1.0-bed_por)*Bed_Frac*(a1*ratio**3+a2*ratio**2&
          &+a3*ratio+a4)
   endif
!JIanzhong

  !Winterwerp
  !Erosion =  MAX(Erate*Bed_Frac*(tau_w/tau_ce-1.0),0.0)
  
  End Function Erosion
!==========================================================================
! Calculate Depositional Flux and Update Concentration from Settling
!==========================================================================
  Subroutine Calc_Deposition
  use all_vars
  implicit none
  integer  :: ised ,i,mcyc,ncyc,k
  real(sp) :: wset(kbm1) ,c(kbm1),dx(kbm1),flux,DTmax,eps,DTdep
  real(sp) :: sal(kbm1),g1(kbm1),ds50,je,tau_cd,cb
  integer  :: cnt,jj,cc
  real(sp),dimension(0:mt,1:kb) :: unode,vnode
  real(sp),dimension(0:mt)      :: cbcn
  real(sp) :: ud


  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Calc_Deposition "

  ! Compute bottom velocity at node
  if(.not.oned_model)then
  CALL E2N3D(u,unode)
  CALL E2N3D(v,vnode)
  endif

  eps = epsilon(eps)

  !Loop over Sediment Class       
  do ised=1,nsed
    do i=1,m
      if(bottom(i,bflag) > 0.0)then
  
      ! no calculation when dry nodes
# if defined (WET_DRY)
      if(iswetn(i)/=1)cycle  
# endif

      !initialize flux
      flux = 0.0

      !setup 1D concentration and grid arrays
      c(1:kbm1) = sed(ised)%conc(i,1:kbm1)
      dx(1:kbm1) = dz(i,1:kbm1)*d(i)
      
      !calculate the settling velocity if sediment is cohesive
      if(sed(ised)%stype == 'cohesive')then
!        call calc_wset(kbm1,c,wset,sed(ised)%Wset)   !Mahta's method
        ds50 = sed(ised)%sd50
        sal(1:kbm1)=s1(i,1:kbm1)
      !calc taub from upwind or downwind elements - beta
# if defined (WAVE_CURRENT_INTERACTION)
      if(MB_BBL_USE.or.SG_BBL_USE.or.SSW_BBL_USE)then
         taub(i) = 1025.0 * taucwmax(i)
      else
# endif
        cnt = 0
        je  = 0.
        taub(i) = 0.
        do jj=1,ntve(i)
           cc=nbve(i,jj)
          !if(xc(cc) < vx(i))then !cell is upstream
          !if(xc(cc) > vx(i))then !cell is downtream
            cnt = cnt + 1
            taub(i) = taub(i) + sqrt(wubot(cc)**2 + wvbot(cc)**2) 
          ! endif
        end do
        taub(i) = 1025.*taub(i)/cnt
# if defined (WAVE_CURRENT_INTERACTION)
        end if
# endif
        g1 = 0.0
        do k=1,kbm1
          if(d(i)<=5.0)then
            je = mf**2*(unode(i,k)**2 + vnode(i,k)**2)/5**(4./3.) 
          else
            je = mf**2*(unode(i,k)**2 + vnode(i,k)**2)/d(i)**(4./3.) 
          end if          
          g1(k)=sqrt(grav*sqrt(unode(i,k)**2+vnode(i,k)**2)*je/Vcons)
        end do
        call calc_wset_floc(kbm1,c,sal,g1,ds50,sed(ised)%w0,wset)
!        if(d(i)<1.0)then
!          ud = 0.812*(ds50*1000.0)**(0.4)*(wset(kbm1)*1000.0*1.0)**(0.2)
!        else
!          ud = 0.812*(ds50*1000.0)**(0.4)*(wset(kbm1)*1000.0*d(i))**(0.2)
!        end if
!        tau_cd = 1025*cbcn(i)*ud**2
        tau_cd = sed(ised)%t_cd(i)
      else
        wset = sed(ised)%Wset
        if(SEDIMENT_PARAMETER_TYPE/=UNIFORM)then
           tau_cd = sed(ised)%t_cd(i)
        else
           tau_cd = sed(ised)%tau_cd
        end if
      endif

      !set up cycles (use max(CFL) == 1)
      DTmax = Settle_CFL*minval(dx(1:kbm1)/wset(1:kbm1))
      mcyc = int(DTsed/DTmax + 1. - eps)
      DTdep = DTsed/float(mcyc)

      !call flux-limited settling equation
      do ncyc = 1,mcyc
        if(sed(ised)%stype == 'cohesive' )then
          call settle_flux_floc(kbm1,c,dx,wset,flux,DTdep,tau_cd,taub(i))
        else
          call settle_flux(kbm1,c,dx,wset,flux,DTdep)
        end if
      end do

! Jianzhong Ge 03/05/2013
      sed(ised)%wset0(i,1:kbm1)=wset(1:kbm1)
! Jianzhong Ge 03/05/2013

      !store solution in 3D array
      sed(ised)%conc(i,1:kbm1) = c(1:kbm1)

      !store depositional flux
      sed(ised)%dflx(i) = flux 

      !shutdown erosion if node is inactive 
      sed(ised)%dflx(i) = sed(ised)%dflx(i)
    else
      sed(ised)%dflx(i) = 0.0
    endif
    end do
  end do

! J. Ge for fluid mud layer
# if defined (FLUID_MUD)
  Settling = 0.0_SP
  do ised=1,nsed
    do i=1,m
# if defined (WET_DRY)
      if(iswetn(i)==1)then
# endif
        Settling(i) = Settling(i) + sed(ised)%dflx(i)
# if defined (WET_DRY)
      else
        Settling(i) = 0.0_sp
      end if
# endif
    end do
  end do  
# endif

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Calc_Deposition "

  End Subroutine Calc_Deposition

!==========================================================================
! Calculate settling velocity of cohesive sediment using the method of
!           Qingyu Sha (including floccuation process)
! Refer : 1.Pingxing Ding et al, 2003, Numerical simulation of total
!         sediment under wave and current in Changjiang Estuary,
!         Acta Oceanologica Sinica, pp114-124,Vol.25, No.5 
!         2.Kelin Hu,2003,2-D Suspended Sediment numerical Simulation
!           under waves and currents in Yangtze Estuary,PHD doctoral
!           dissertation 2003 East China Normail University
!    or
!    User Defined
!==========================================================================
  Subroutine Calc_Wset_Floc(n,c,s,g1,ds50,w0,wset) 
  implicit none
  integer , intent(in   ) :: n
  real(sp), intent(in   ) :: ds50,w0, c(n),s(n),g1(n)
  real(sp), intent(inout) :: wset(n)
 
  integer :: i

  !setup for Changjiang Estuary
  real(sp) :: a = 0.274 ! proportional constant
  real(sp) :: b = 0.48  ! exponential constant of suspended sediment concertration
  real(sp) :: d = 0.03  ! exponential constant of salinity
  real(sp) :: e = 0.22  ! exponential constant of turblence intensity
  real(sp) :: f = 0.58  ! exponential constant of Ds50
  real(sp) :: factor

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Calc_Wset_Floc"

  do i=1,n
     factor = a*c(i)**b*s(i)**d*g1(i)**e/((ds50*1000.0)**f)
     if(factor<1.0.or.s(i)<5.0.or.s(i)>15.0)factor = 1.0
     wset(i)=w0*factor
     if(wset(i)>0.5e-3)wset(i)=0.5e-3
  end do

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Calc_Wset_Floc"

  End Subroutine Calc_Wset_Floc

  Subroutine Calc_Wset_hindered(n,c,w0,wset,cmax,scale) 
  implicit none
  integer , intent(in   ) :: n
  real(sp), intent(in   ) :: w0,cmax,scale,c(n)
  real(sp), intent(inout) :: wset(n)
 
  integer :: i

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Calc_Wset_Floc2"

  do i=1,n
    if(c(i) < cmax)then
       wset(i) = w0
    else
      wset(i)=w0*(1.0-c(i)/1650.0)**5
      !wset(i) = w0*scale
      !wset(i) = w0*(cmax/c(i))**3
    endif
  end do


  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Calc_Wset_Floc2"

  End Subroutine Calc_Wset_hindered

!Jianzhong Ge 03/05/2013
  Subroutine Hindered_Diffusion
  use lims, only: m
  use all_vars,only: kh,kbm1
  implicit none
  real(sp) :: sed_conc(1:m,1:kbm1)  !mass concentration
  real(sp) :: vol_conc(1:m,1:kbm1)  !volume concentration
  integer :: i,ii,k
  real,parameter :: cs0=0.65
  real :: phi

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Hindered_Diffusion"

  sed_conc = 0.0  
  do ii=1,nsed
     sed_conc=sed_conc+sed(ii)%cnew
  end do

  do i=1,m
    do k=1,kbm1
      if(sed_conc(i,k)>2.5)then
!        vol_conc(i,k)=sed_conc(i,k)/2650.0
!        phi=1+(vol_conc(i,k)/Cs0)**0.8-2.0*(vol_conc(i,k)/Cs0)**0.4
!        kh(i,k)=kh(i,k)*phi**4.0
         kh(i,k)=kh(i,k)*0.02
!        print*,kh(i,k),phi**4.0,vol_conc(i,k),sed_conc(i,k)
      end if
    end do
  end do

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Hindered_Diffusion"

  End Subroutine Hindered_Diffusion
!Jianzhong Ge 03/05/2013

!==========================================================================
! Calculate settling velocity of cohesive sediment using the method of
!    Mehta, 1986 or
!    Hindered settlement strategy
!    or
!    User Defined
!==========================================================================
  Subroutine Calc_Wset(n,c,wset,Wset_Mean) 
  implicit none
  integer , intent(in   ) :: n
  real(sp), intent(in   ) ::  c(n)
  real(sp), intent(inout) :: wset(n)
  real(sp), intent(in   ) :: Wset_Mean
  integer :: i

  !setup for Ariake Bay
  real(sp) :: a = 0.5 !Mehta proportional constant
  real(sp) :: b = 2.0 !Mehta exponential constant
  real(sp) :: cmax = 1.0 !Concentration for hindered settling

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Calc_Wset"

  do i=1,n
    if(c(i) < cmax)then
       wset(i) = a*(c(i)**b)  
    else
       !hindered settling here
    endif
  end do

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Calc_Wset"

  End Subroutine Calc_Wset




!==========================================================================
! Calculate settling velocity of cohesive sediment using the method of
!    Richardson & Zaki (1954)
!    Hindered settlement strategy
!    or
!    User Defined
!==========================================================================
  Subroutine Calc_Wset_delft(n,c,Wset_max,wset) 
  implicit none
  integer , intent(in   ) :: n
  real(sp), intent(in   ) :: c(n)
  real(sp), intent(inout) :: wset(n)
  real(sp), intent(in   ) :: Wset_Max
  integer :: i

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Calc_Wset_delft"

  do i=1,n
     wset(i)=(1-c(i)/1600.0)**5.0*Wset_max
  end do

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Calc_Wset_delft"

  End Subroutine Calc_Wset_delft

!==========================================================================
! Calculate Settling Flux and Update Sediment Concentration in Water Column 
!   use the SLIP limiter
!   Jameson, A., Analysis and Design of Numerical Schemes for Gas Dynamics I
!      Artificial Diffusion, Upwind Biasing, Limiters, and their effects
!      on Accuracy and Multigrid Convergence, International Journal of 
!      Computational Fluid Dynamics, Vol 4, 171-218, 1995.
!
!      -second order accurate in smooth regions
!      -guaranteed local extremum diminishing
!      
!   conv:  convective flux
!   diss:  dissipative flux
!   
!==========================================================================
  Subroutine Settle_Flux(n,c,dx,wset,flux,deltat) 
  implicit none
  integer , intent(in   ) :: n
	  real(sp), intent(inout) ::  c(n)
	  real(sp), intent(in   ) :: dx(n)
	  real(sp), intent(in   ) :: wset(n)
	  real(sp), intent(inout) :: flux 
	  real(sp), intent(in   ) :: deltat 
	  real(sp) :: conv(n+1),diss(n+1)
	  real(sp) :: cin(-1:n+2)
	  real(sp) :: win(-1:n+2)
	  real(sp) :: dis4,wvel
  integer  :: i

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Settle_Flux "

  !transfer to working array
  cin(1:n) = c(1:n)
  win(1:n) = wset(1:n) 

  !surface bcs (no flux)
  cin(0)  =  -cin(1) 
  cin(-1) =  -cin(1)
  win(0)  =  -wset(1)
  win(-1) =  -wset(1)

  !bottom bcs (extrapolate)
  cin(n+1) = cin(n) 
  cin(n+2) = cin(n)
  win(n+1) = win(n)
  win(n+2) = win(n)

  !flux computation
  do i=1,n+1
    wvel    = .5*(win(i)+win(i-1))  !settle velocity at interface
    dis4    = wvel/2.
    conv(i) = wvel*(cin(i)+cin(i-1))/2. 
    diss(i) = dis4*(cin(i)-cin(i-1)-lim(cin(i+1)-cin(i),cin(i-1)-cin(i-2))) 
  end do

  !zero out surface flux
  conv(1) = 0. ; diss(1) = 0.

  !update
  do i=1,n
    c(i) = cin(i) + (deltat/dx(i))*(-conv(i+1)+conv(i) + diss(i+1)-diss(i)) 
  end do

  !set bottom flux (> 0 = deposition)
  flux = flux + deltat*(conv(n+1)-diss(n+1))

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Settle_Flux "

 End Subroutine Settle_Flux


!==========================================================================
! Calculate Settling Flux and Update Sediment Concentration in Water Column 
!           use cohesive sediment method with floccuation process
! Refer : 1.Pingxing Ding et al, 2003, Numerical simulation of total
!         sediment under wave and current in Changjiang Estuary,
!         Acta Oceanologica Sinica, pp114-124,Vol.25, No.5 
!         2.Kelin Hu,2003,2-D Suspended Sediment numerical Simulation
!           under waves and currents in Yangtze Estuary,PHD doctoral
!           dissertation 2003 East China Normail University      
!
!   conv:  convective flux
!   diss:  dissipative flux
!   
!==========================================================================
  Subroutine Settle_Flux_Floc(n,c,dx,wset,flux,deltat,tau_cd,taub) 
  implicit none
  integer , intent(in   ) :: n
	  real(sp), intent(inout) ::  c(n)
	  real(sp), intent(in   ) :: dx(n)
	  real(sp), intent(in   ) :: wset(n)
	  real(sp), intent(inout) :: flux 
	  real(sp), intent(in   ) :: deltat 
          real(sp), intent(in   ) :: tau_cd,taub
	  real(sp) :: conv(n+1),diss(n+1)
	  real(sp) :: cin(-1:n+2)
	  real(sp) :: win(-1:n+2)
          real(sp) :: wset2(n)
	  real(sp) :: dis4,wvel
 
  integer  :: i

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Settle_Flux_Floc "
  
  !apply the critical shear stress for deposition
  wset2 = wset
!  if(taub>=tau_cd)then
!    wset2(n)=0.0
!  end if

  !transfer to working array
  cin(1:n) = c(1:n)
  win(1:n) = wset2(1:n) 

  !surface bcs (no flux)
  cin(0)  =  -cin(1) 
  cin(-1) =  -cin(1)
  win(0)  =  -wset2(1)
  win(-1) =  -wset2(1)

  !bottom bcs (extrapolate)
  cin(n+1) = cin(n) 
  cin(n+2) = cin(n)
  win(n+1) = win(n)
  win(n+2) = win(n)
  !flux computation
  do i=1,n+1
    wvel    = .5*(win(i)+win(i-1))  !settle velocity at interface
    dis4    = wvel/2.
    conv(i) = wvel*(cin(i)+cin(i-1))/2. 
    diss(i) = dis4*(cin(i)-cin(i-1)-lim(cin(i+1)-cin(i),cin(i-1)-cin(i-2))) 
  end do
  
  !zero out surface flux
  conv(1) = 0. ; diss(1) = 0.
  
  !update
  do i=1,n
    c(i) = cin(i) + (deltat/dx(i))*(-conv(i+1)+conv(i) + diss(i+1)-diss(i)) 
  end do

  !set bottom flux (> 0 = deposition)
  flux = flux + deltat*(conv(n+1)-diss(n+1))

  if(flux<0)call fatal_error("Error: Deposition Flux is negative!!!")

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Settle_Flux_Floc "

 End Subroutine Settle_Flux_Floc





!==========================================================================
! Calculate LED Limiter L(u,v)  
!==========================================================================
  Function Lim(a,b)
  real(sp) lim,a,b
  real(sp) q,R
  real(sp) eps
  eps = epsilon(eps)
  
 ! exponent
 ! qlim = 0. !1st order
 ! qlim = 1. !minmod
 ! qlim = 2. !van leer

  R = abs(   (a-b)/(abs(a)+abs(b)+eps) )**qlim
  lim = .5*(1-R)*(a+b)

  End Function Lim

!==========================================================================
! Update Surface Bed Properties Due to Bedload/SusLoad Flux 
! Change:
!   sed(:)%mass(:,1)   to reflect increased/decrease mass (kg)
!   bed(:,1,itchk)     to reflect changed bed thickness properties
!==========================================================================

  Subroutine Update_Surface_Bed(fluxtype)
  use lims, only: m
  implicit none
  character(len=*) :: fluxtype
  real(sp) :: dz,accum,flux,t_mass,eps,junk
  integer  :: i,ised

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Update_Surface_Bed "

  eps = epsilon(eps)

  !======= adjust for bedload flux ===================
  if(fluxtype == 'bedload')then

  !top layer mass and thickness
  !bflx > 0 means material is lost from the bed.
  do i=1,m
    do ised=1,nsed
      flux = sed(ised)%bflx(i)

      !change top layer bed mass of each sed type 
      accum  = sed(ised)%mass(i,1) - flux 
      sed(ised)%mass(i,1) = max(accum,0.0) 

      !change top layer bed thickness 
      dz = flux/(sed(ised)%srho*(1.0-bed(i,1,iporo)))
      junk = bed(i,1,ithck)
      bed(i,1,ithck) = max(bed(i,1,ithck)-dz,0.0)
      bottom(i,morph) = bottom(i,morph) - dz 

    end do
  end do

  !top layer fractions 
  do i=1,m
    t_mass = 0.0_sp
    do ised=1,nsed
      t_mass = t_mass + sed(ised)%mass(i,1)  
    end do
    do ised=1,nsed
      sed(ised)%frac(i,1) = sed(ised)%mass(i,1)/MAX(t_mass,eps)
    end do
  end do

  !======= adjust for susload flux ===================
  elseif(fluxtype =='susload')then
  
  do i=1,m
    do ised=1,nsed
      flux = (sed(ised)%eflx(i)-sed(ised)%dflx(i))*morpho_factor
      !if depositional + first deposit time step, save mass and loop 
      !http://woodshole.er.usgs.gov/project-pages/sediment-transport/ 
      sed(ised)%depm(i) = 0.0
      if(flux < 0)then                                  !!depositional
        if(T_model > bed(i,1,iaged)+1.1*DTsed .and. &   !!first deposit time step
           bed(i,1,ithck) > thck_cr)then                !!thickness surpasses critical
          sed(ised)%depm(i) = -flux
          cycle
        else                                            !!update age of surface layer
          bed(i,1,iaged) = T_model
        endif
      endif

      !change top layer bed mass of each sed type 
      accum  = sed(ised)%mass(i,1) - flux 
      sed(ised)%mass(i,1) = max(accum,0.0) 

      !change top layer bed thickness 
      dz = flux/(sed(ised)%srho*(1.0-bed(i,1,iporo)))
      bed(i,1,ithck) = max(bed(i,1,ithck)-dz,0.0)
      bottom(i,morph) = bottom(i,morph) - dz 

    end do
  end do

  !top layer fractions 
  do i=1,m
    t_mass = 0.0_sp
    do ised=1,nsed
      t_mass = t_mass + sed(ised)%mass(i,1)  
    end do
    do ised=1,nsed
      sed(ised)%frac(i,1) = sed(ised)%mass(i,1)/MAX(t_mass,eps)
    end do
  end do

  !======= error ======================
  else
    write(*,*)'argument to Update_Surface_Bed must be: "susload" or "bedload"'
    stop
  endif !fluxtype

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Update_Surface_Bed "

  End Subroutine Update_Surface_Bed


!==========================================================================
! Update Bottom Properties due to changes in surface layer sed fractions 
! Use geometric mean, good for log normal distribution
!==========================================================================
  Subroutine Update_Bottom_Properties
  use lims, only: m
  implicit none
  integer  :: i,ised
  real(sp) :: temp,sum1,sum2,sum3,sum4,eps

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Update_Bottom_Properties "

  ! initialize
  eps = epsilon(temp) 
  
  ! update bottom properties using new bed fractions
  do i=1,m
    sum1 = 1.0 ; sum2 = 1.0 ; sum3 = 1.0 ; sum4 = 1.0
    do ised=1,nsed
      if(SEDIMENT_PARAMETER_TYPE/=UNIFORM)then
        sum1 = sum1*(sed(ised)%t_ce(i))**sed(ised)%frac(i,1) 
      else
        sum1 = sum1*(sed(ised)%tau_ce)**sed(ised)%frac(i,1)
      end if  
      sum2 = sum2*(sed(ised)%Sd50   )**sed(ised)%frac(i,1)
      sum3 = sum3*(sed(ised)%wset   )**sed(ised)%frac(i,1)
      sum4 = sum4*(sed(ised)%Srho   )**sed(ised)%frac(i,1)
    end do
    bottom(i,itauc) = sum1
    bottom(i,isd50) = sum2
    bottom(i,isphi) = -log(bottom(i,isd50)*1000.)/log(2.0)
    bottom(i,iwset) = sum3
    bottom(i,idens) = max(sum4,min_Srho) 
  end do 
  
  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Update_Bottom_Properties "

  End Subroutine Update_Bottom_Properties 

!==========================================================================
! Update Active Layer Thickness 
!   Use method of Harris and Wiberg (1997):
!   -Approaches to quantifying long-term continental shelf sediment transport
!   -with an example from the northern California STRESS mid-shelf site.
!   [Continental Shelf REsearch, 17, 1389-1418]
!
!   z_a = max[k_1*(tau_bot - tau_ce)*rho_0 , 0] + k_2 * D_50
!       where
!       k_1 = .007 (empirical constant)
!       k_2 = 6.0  (empirical constant)
!       tau_bot = shear stress on the bed
!       tau_ce  = critical shear stress for erosion of the bed (averaged over
!                 sediment classes)
!       D_50    = median grain diameter of the surface sediment
!       note: our shear stresses are nondimensionalized by rho, so
!             rho does not appear in actual calculation (below)
!==========================================================================
  Subroutine Calc_Active_Layer 
  use lims,     only: m
  implicit none
  integer  :: i
  real(sp) :: mfac

  mfac = dim(morpho_factor,1.0) 
  !calculate active layer thickness 
  do i=1,m
    bottom(i,iactv) = mfac*max(0.,.007*(taub(i)-bottom(i,itauc)))+6.0*bottom(i,isd50)
  end do

  End Subroutine Calc_Active_Layer
  
!==========================================================================
! Update Stratigraphy       
!   -Calculate empirical active layer thickness from bottom stress
!   -Expand top layer to this thickness by shifting necessary mass up
!    through the bed.
!==========================================================================

  Subroutine Update_Stratigraphy
  use lims,     only: m
  implicit none
  integer :: i,k,ksed,ks,ised
  real(sp) :: thck_avail,thck_to_add,eps,tmp1,tmp2,tmp3,top_layer_mass
  real(sp) :: delta(m) 
 
  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Update_Stratigraphy "

  !initialize
  eps = epsilon(tmp1) 
  delta = 0.0

  !calculate deficit between active layer thickness and surface layer thickness
  do i=1,m
    delta(i) = bottom(i,iactv)-bed(i,1,ithck)
  end do
  
  ! (single layer system) ensure top layer > active layer thickness 
  if(nbed == 1)then 
    do i=1,m
       if(delta(i) > 0) bottom(i,iactv) = bed(i,1,ithck)
    end do
    return
  endif
   
  ! (multi layer system) ensure top layer > active layer thickness 
  do i=1,m
    thck_to_add = delta(i)
    if(thck_to_add > 0.0)then !must redistribute layers 
      thck_avail=0.0

      !find fractional layer (below which there will be no mass after forming active layer) 
      Ksed=1 
      do k=2,nbed
        if (thck_avail < thck_to_add) then
          thck_avail=thck_avail+bed(i,k,ithck)
          Ksed=k
        end if
      end do

      !if not enough material in bed to satisfy active layer req, use all available 
      if (thck_avail < thck_to_add) then
        bottom(i,iactv)=bed(i,1,ithck)+thck_avail
        thck_to_add=thck_avail
      end if
  
      !update the bed mass of surface and fractional layers
      tmp2=MAX(thck_avail-thck_to_add,0.0)/MAX(bed(i,Ksed,ithck),eps)
      do ised=1,nsed
        tmp1=0.0 
        do k=1,Ksed
          tmp1=tmp1+sed(ised)%mass(i,k)
        end do
        sed(ised)%mass(i,1   ) = tmp1-sed(ised)%mass(i,Ksed)*tmp2
        sed(ised)%mass(i,Ksed) = sed(ised)%mass(i,Ksed)*tmp2
      end do

      !update thickness of fractional layer Ksed
      bed(i,Ksed,ithck)=MAX(thck_avail-thck_to_add,0.0)

      !update top layer bed fraction 
      top_layer_mass = 0.0
      do ised=1,nsed
        top_layer_mass = top_layer_mass + sed(ised)%mass(i,1)
      end do
      do ised=1,nsed
        sed(ised)%frac(i,1)=sed(ised)%mass(i,1)/MAX(top_layer_mass,eps) 
      end do

      !update bed thickness of top layer
      bed(i,1,ithck)=bottom(i,iactv)

      !pull layers Ksed to Bottom up to fill layer 2 down 
      do k=Ksed,Nbed
        ks=Ksed-2
        bed(i,k-ks,ithck)=bed(i,k,ithck)
        bed(i,k-ks,iporo)=bed(i,k,iporo)
        bed(i,k-ks,iaged)=bed(i,k,iaged)
        do ised=1,nsed
          sed(ised)%frac(i,k-ks) = sed(ised)%frac(i,k)
          sed(ised)%mass(i,k-ks) = sed(ised)%mass(i,k)
        end do
      end do
 
      !split what was in the bottom layer to fill empty bottom cells
      !note:  porosity of nbed (bed(i,nbed,iporo)) does not change.
      ks=Ksed-2
      tmp3=1.0/real(ks+1)
      do k=Nbed,Nbed-ks,-1
        bed(i,k,ithck)=bed(i,Nbed-ks,ithck)*tmp3
        bed(i,k,iaged)=bed(i,Nbed-ks,iaged)
        do ised=1,nsed
          sed(ised)%frac(i,k)=sed(ised)%frac(i,Nbed-ks)
          sed(ised)%mass(i,k)=sed(ised)%mass(i,Nbed-ks)*tmp3
        end do
      end do

     end if !thck_to_add > 0
  end do !loop over nodes

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Update_Stratigraphy "
    
  End Subroutine Update_Stratigraphy

!==========================================================================
! Add New Layer                 
!   If first deposit time step:  
!       a.) combine bottom layers 
!       b.) shift upper layers down
!       c.) create new top layer
!==========================================================================

  Subroutine Add_New_Layer 
  use lims,     only: m
  implicit none
  integer :: i,k,ksed,ised
  real(sp) :: eps,t_mass,cnt
 
  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Add_New_Layer "

  !initialize
  eps = epsilon(eps) 

  do i=1,m
  
    cnt = 0.0
    do ised=1,nsed
      cnt = cnt + sed(ised)%depm(i)
    end do

    if(cnt > 0)then !!need new surface layer

      if(nbed > 1)then!--->
      !combine bottom two layers
      bed(i,nbed,ithck) =     bed(i,nbed-1,ithck) + bed(i,nbed,ithck)
      bed(i,nbed,iporo) = .5*(bed(i,nbed-1,iporo) + bed(i,nbed,iporo))
      bed(i,nbed,iaged) = .5*(bed(i,nbed-1,iaged) + bed(i,nbed,iaged))

      t_mass =0.0
      do ised=1,nsed
        sed(ised)%mass(i,nbed) = sed(ised)%mass(i,nbed-1)+sed(ised)%mass(i,nbed) 
        t_mass = t_mass +sed(ised)%mass(i,nbed)
      end do
      do ised=1,nsed
        sed(ised)%frac(i,nbed) = sed(ised)%mass(i,nbed)/MAX(t_mass,eps)
      end do

      !push layers down
      do k=Nbed-1,2,-1
        bed(i,k,ithck) = bed(i,k-1,ithck)
        bed(i,k,iporo) = bed(i,k-1,iporo)
        bed(i,k,iaged) = bed(i,k-1,iaged)
        do ised =1,nsed
          sed(ised)%frac(i,k) = sed(ised)%frac(i,k-1)
          sed(ised)%mass(i,k) = sed(ised)%mass(i,k-1)
        end do
      end do

      !refresh top layer in multi layer bed
      bed(i,1,ithck)=0.0
      do ised=1,nsed
        sed(ised)%mass(i,1)=0.0
      end do
      end if ! <---

     

      !update surface layer properties
      t_mass = 0.0
      do ised=1,nsed
        sed(ised)%mass(i,1) = sed(ised)%mass(i,1)+ sed(ised)%depm(i)
        t_mass = t_mass + sed(ised)%mass(i,1)
        bed(i,1,ithck)=bed(i,1,ithck)+sed(ised)%depm(i)/    &
                         (sed(ised)%Srho*(1.0-bed(i,1,iporo)))
      end do
      bed(i,1,iaged) = T_model  
      do ised=1,nsed
        sed(ised)%frac(i,1) = sed(ised)%mass(i,1)/MAX(t_mass,eps) 
        sed(ised)%depm(i  ) = 0.0
      end do
 
    end if  !need new surface layer
  end do !node loop

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Add_New_Layer "

  End Subroutine Add_New_Layer

!==========================================================================
! Update Total Sediment Thickness for Reporting
!==========================================================================

  Subroutine Update_Thickness_Delta
  use lims,     only: m
  implicit none
  integer :: i,k,ised
 
  !shift last thickness to lthck
  bottom(:,lthck) = bottom(:,nthck)
  bottom(:,nthck) = 0.0

  !calculate new thickness
  do i=1,m
    do k=1,nbed
      bottom(i,nthck) = bottom(i,nthck) + bed(i,k,ithck)
    end do
  end do

  !calculate new delta
  do i=1,m
    bottom(i,dthck) = bottom(i,dthck) + bottom(i,nthck) - bottom(i,lthck)
  end do

  !calculate new total mass
  do i=1,m
    bottom(i,tmass) = 0.0 
    do ised=1,nsed
      do k=1,nbed
        bottom(i,tmass) = bottom(i,tmass) + sed(ised)%mass(i,k)
      end do
    end do
  end do
  
  End Subroutine Update_Thickness_Delta

!==========================================================================
! Report Layer Stats for Debug   
!==========================================================================

  Subroutine Layer_Report(i)
  implicit none
  integer, intent(in) :: i
  integer :: j,ised

  write(*,*)'----------bottom props'
  write(*,*)'active thickness',bottom(i,iactv)
  write(*,*)'total  thickness',bottom(i,nthck)
  write(*,*)'mean grain size ',bottom(i,isd50)
  write(*,*)'mean density    ',bottom(i,idens)
  write(*,*)'mean settle V   ',bottom(i,iwset)
  write(*,*)'mean E stress   ',bottom(i,itauc)

  write(*,*)'--------- bed props'
  write(*,*)'total  thickness',bottom(i,nthck)
  write(*,*)'total  mass'     ,bottom(i,tmass)
  if(bottom(i,nthck) <= 0.)stop

  End Subroutine Layer_Report

!==========================================================================
!  Setup open boundary nudging for sediment obc
!   - check for existence of sed nudging file
!   - read data (concentration of each sed type at each obc)
!==========================================================================
  Subroutine Setup_Sed_OBC
  Use Control,  only : casename,input_dir,msr,serial,par
!  Use Mod_OBCs, only : iobcn,i_obc_gl,iobcn_gl
  Use Mod_Utils, only: pstop
# if defined (MULTIPROCESSOR)
  Use Mod_Par,  only : nlid
# endif

  integer :: i,j,tmp,ncnt,i1
  integer,  allocatable :: vtmp(:),atemp(:,:)
  real(sp), allocatable :: sed_obc_gl(:,:)
  integer, parameter    :: sedobc = 81
  character(len=120)    :: fname
  logical :: fexist
  

  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Setup_Sed_OBC "

  !------------------------------------------------------
  ! make sure sediment nudging data file exists 
  !------------------------------------------------------
  fname = "./"//trim(input_dir)//"/"//trim(casename)//"_sed_nudge.dat"

  inquire(file=trim(fname),exist=fexist)
  if(msr .and. .not.fexist)then
    write(*,*)'sediment nudging data: ',fname,' does not exist'
    write(*,*)'halting.....'
    call pstop
  end if

  open(unit=sedobc,file=fname,form='formatted')
  !------------------------------------------------------
  ! allocate global data
  !------------------------------------------------------
  allocate(sed_obc_gl(Nobc_gl,nsed)) ; sed_obc_gl = 0.

  !header
  read(sedobc,*)
  read(sedobc,*)
  read(sedobc,*)

  !# obc nodes
  read(sedobc,*) tmp 

  if(tmp /= Nobc_gl .and. msr)then
    write(*,*)'number of open boundary nodes for sediment nudging: ',tmp
    write(*,*)'this is not consistent with number of open boundary nodes in'
    write(*,*)'the model: ',Nobc_gl
    write(*,*)'halting'
    call pstop
  endif

  allocate(vtmp(Nobc_gl))
  do i=1,Nobc_gl  
    read(sedobc,*) i1,vtmp(i),(sed_obc_gl(i,j),j=1,nsed)
  end do
  close(sedobc)

  ! !------------------------------------------------------
  ! ! make sure global obc node list is consistent
  ! !------------------------------------------------------
  ! 
  ! if(msr)then
  ! do i=1,iobcn_gl
  !   if(i_obc_gl(i) /= vtmp(i))then
  !     write(*,*)'obc node list for sediment nudging not consistent with model'
  !     write(*,*)'obc node list' 
  !     write(*,*)'halting'
  !     call pstop
  !   end if
  ! end do
  ! endif
    
  !---------------------------------------------------------
  ! shift sediment nudging data to local open boundary nodes
  !---------------------------------------------------------

  if(serial)then
    do i=1,Nobc_gl
      do j=1,nsed
        sed(j)%cobc(i) = sed_obc_gl(i,j)
      end do
    end do
  endif

# if defined (MULTIPROCESSOR)
  if(par)then
    allocate(atemp(Nobc_gl,nsed))
    ncnt = 0
    !!set up local open boundary nodes
    do i=1,Nobc_gl
      i1 = nlid( vtmp(i) )
      if(i1 /= 0)then
        ncnt = ncnt + 1
        do j=1,nsed
          atemp(ncnt,j) = sed_obc_gl(i,j)
        end do
      end if
    end do

    !transfer sed nudging data global --> local
    if(ncnt > 0)then
      do j=1,nsed
        sed(j)%cobc(1:ncnt) = atemp(1:ncnt,j)
      end do
    end if

    deallocate(atemp) 
  end if
# endif

  if(msr)write(*,*)'obc setup for sediment model: complete'

  !---------------------------------------------------------
  ! cleanup 
  !---------------------------------------------------------
  deallocate(sed_obc_gl)

  if(dbg_set(dbg_sbr)) write(ipt,*) "End: Setup_Sed_OBC "

  End Subroutine Setup_Sed_OBC



   SUBROUTINE UPDATE_DEPTH
!------------------------------------------------------------------------------|
   USE ALL_VARS
   USE MOD_OBCS
   IMPLICIT NONE
   REAL(SP) :: TEMP, DX, DY
   INTEGER  :: I,I1,K,J1,J2,CC
   LOGICAL  :: ISONB2
!------------------------------------------------------------------------------|
   IF(DBG_SET(DBG_SBR)) write(ipt,*) "START: UPDATE_DEPTH"

# if defined (MULTIPROCESSOR)
    !IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,bottom )
    !IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,H )
# endif

   DO I = 1,M
!
!   DEPTHS ON OPEN BOUNADY NODES AND THE NODES NEXT TO THOSE NODES ARE KEPT INVARIANT
!
    ISONB2=.FALSE.
     DO I1 = 1,IOBCN
       J1 = I_OBC_N(I1)
       J2 = NEXT_OBC(I1)
       IF(I==J1 .OR. I==J2) ISONB2=.TRUE.
     ENDDO
    IF(.NOT. ISONB2) H(I) = H(I) - bottom(I,morph)
    bottom(I,morph) = 0.0
   ENDDO
# if defined (MULTIPROCESSOR)
    IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,H )
# endif
   CALL N2E2D(H,H1)
   !DO I = 1, N
   ! CC = 0
   ! H1(I) = 0.
   ! DO K = 1,3
   !  DX = VX(NV(I,K)) - XC(I)
   !  DY = VY(NV(I,K)) - YC(I)
   !  IF( (DX*U(I,KBM1) - DY*V(I,KBM1)) <  0.) THEN !node is upstream
   !   CC = CC + 1
   !   H1(I) = H1(I) + H(NV(I,K))
   !  ENDIF
   ! ENDDO 
   !  H1(I) = H1(I)/CC
   !ENDDO
# if defined (MULTIPROCESSOR)
    !IF(PAR)CALL AEXCHANGE(EC,MYID,NPROCS,H1 )
    !IF(PAR)CALL AEXCHANGE(EC,MYID,NPROCS,EL1,ET1)
    !IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,EL ,ET )
# endif
   D   = H + EL
   DT  = H + ET
   D1  = H1+EL1
   DT1 = H1+ET1
# if defined (MULTIPROCESSOR)
    !IF(PAR)CALL AEXCHANGE(EC,MYID,NPROCS,D1,DT1 )
    !IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,D ,DT )
# endif


   IF(DBG_SET(DBG_SBR)) write(ipt,*) "END: UPDATE_DEPTH"

   END SUBROUTINE UPDATE_DEPTH

# endif
   

End Module Mod_Sed 
